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Abstract 

The equation of state of QCD at finite temperatures and baryon densities has a wide range 
of apphcations in many fields of modern particle and nuclear physics. It is the main ingredient 
to describe the dynamics of experimental heavy ion collisions, the expansion of the early universe 
in the standard model era and the interior of compact stars. On most scales of interest, QCD 
is strongly coupled and not amenable to perturbative investigations. Over the last decade, first 
principles calculations using lattice QCD have reached maturity, in the sense that for particular 
discretisation schemes simulations at the physical point have become possible, finite temperature 
results near the continuum limit are available and systematic errors begin to be controlled. This 
review summarises the current theoretical and numerical state of the art based on staggered and 
Wilson fermions. 
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1 Introduction 



Quantum Chromodynamics (QCD) is the fundamental quantum field theory describing the strong 
interactions within the Standard Model. As such it is also the fundamental theory of nuclear matter. 
One of its many well-tested features is asymptotic freedom, according to which the coupling vanishes 
at asymptotically high energy scales revealing the nature of its constituents, the quarks and gluons. 
Conversely, the coupling is strong on hadronic energy scales < 1 GeV and we observe confinement of 
quarks and gluons which is not amenable to weak coupling expansions. Asymptotic freedom implies 
fascinating changes of dynamics when QCD is considered under extreme conditions, where either a high 
temperature sets the dominant momentum scale, as in the early universe, or a high baryon density, 
as in compact stars. Heavy ion collision experiments have succeeded in creating the hot quark gluon 
plasma in the laboratory and operate to further understand its properties, while future heavy ion 
experiments and astronomical observations aim to investigate cold and dense matter. Of particular 
phenomenological importance is the equation of state of the hot and/or dense system. As a fundamental 
property of the thermal system, it is in principle accessible to experiment allowing for a direct comparison 
with theoretical predictions. Moreover, the equation of state constitutes important input for further 
dynamical analyses of thermal systems, such as the hydrodynamical description of heavy ion collisions [T] 
or the search for new physics in the early universe [2]. 

Because of the inherently non-perturbative nature of the theory, numerical simulations of lattice 
QCD are the only known tool allowing for predictions from first principles, for which it is known 
how to remove the associated systematic errors. Investigations of the equation of state have been 
going on for about two decades. After rapid initial successes with the pure gauge plasma, the step to 
include dynamical fermions has proved soberingly difficult. Systematic errors associated with fermion 
discretisations were initially underestimated, leading to apparent contradictions. However, after an 
impressive collective effort in man and machine power, these issues appear to be finally resolved. It 
was demonstrated that, while it may take some time and effort, systematic errors eventually can be 
controlled and removed, teaching us a lot about the underlying dynamics in the process. The equation 
of state at finite temperature and zero baryon density is now known for Nj = 2 + 1 quark fiavours with 
physical masses very close to its continuum limit, which should be taken within the next year. This is 
based on the staggered fermion discretisation. Calculations with Wilson fermions are somewhat behind, 
but will soon serve as an independent cross check for remaining theoretical issues with the staggered 
formulation. Refinements to include the charm quark are also well on the way. The situation at finite 
density is much less satisfactory. Because of the sign problem of lattice QCD, direct simulations at 
finite baryon density are impossible and further approximations have to be made which are valid for 
sufficiently small chemical potentials. Nevertheless, important first steps in this direction have been 
made and we also understand the response of the equation of state to a small baryon chemical potential. 

Rather than on a complete history, the focus in this review is on the more recent calculations closest 
to the continuum or the physical point, and hence on staggered and Wilson fermions. The control 
of systematic errors being the main task of current lattice calculations. Sees. [3l|5] are devoted to a 
detailed discussion of the different discretisation schemes and their associated cut-off effects along with 
improvement schemes to minimise them. The aim is to also give the non-practitioner some insight into 
what is being done in different calculations. The numerical results on the equation of state are then 
collected in SecsJMTUl 

2 Statistical mechanics of QCD 

We wish to describe a system of particles in some volume V which is in thermal contact with a heat 
bath at temperature T. Associated with the particles may be a set of conserved charges iVj, i = 1, 2, . . . 



3 



(such as particle number, electric charge, baryon number etc.). In quantum field theory, the most direct 
description is in terms of the grand canonical ensemble. Its density operator and partition function are 
given as 

p^Q-W-t^i^i) ^ Z^trp, tr(...) = J](n|(...)|n) , (1) 

n 

where are chemical potentials for the conserved charges, and the quantum mechanical trace is a sum 
over all energy eigenstatcs of the Hamiltonian. Thermodynamic averages for an observable O are then 
obtained as (O) = Z~^Tr(pO) . 

From the partition function, all other thermodynamic equilibrium quantities follow by taking appro- 
priate derivatives. In particular, the (Helmholtz) free energy, pressure, entropy, mean values of charges 
and energy are obtained as 

F = -TlnZ, d{T\nZ) 
a(TlnZ) = ' 

a(TlnZ)' ^ - -pvItS + ,.N.. (2) 



S 



dT 



Since the free energy is known to be an extensive quantity, F — fV, and we are interested in the 
thermodynamic limit, it is often more convenient to consider the corresponding densities, 

f F S N, E 

^^y' '^y' ""'^v' '^y' 

for which we in addition have 

- 'F'^ 9\nZ _ 1 dlnZ 

In a relativistic system, another quantity of interest is the energy momentum tensor for a liquid or gas, 

T'^'^ = (p + e)u''u'' - pg''" , (5) 

where u'^ = (1,0,0,0) is the four velocity of a volume element in the rest frame of the medium. Its 
so-called trace anomaly corresponds to 

I{T) ^ T^^iT) = T'-^P^ . (6) 

Various useful relations between these quantities are then 

TO e + p 2 dp 

7 = e - 3p , s = , = — , (7) 

where the energy derivative of the pressure in the last equation corresponds to the speed of sound in 
the thermal medium. The equation of state is the functional relationship among the thermodynamic 
parameters for a system in thermal equilibrium, 

f{p,V,T,f,,)^0, (8) 

and constitutes the input for many further analyses of a thermal system. 



4 



2.1 QCD at finite temperature and quark density 

Let us now consider the grand canonical partition function of QCD. The derivation of its path integral 
representation is discussed in detail in the textbooks [3J, 

Z{V, i^f, T- g, ruf) = tr (e-(^-^/«/)/^) = j DAD^j Dijj e-^^'^^l ^-Sf[^,^,A,]^ (g) 

with the Euclidean gauge and fermion actions 

l/T 

Sg[A^] = J dr J £x^TtF^,{x)F^,{x), 

V 

Sf[i!,tp,A^] = j dr j d^x^i}f{x){'^^D^ + mf - iipQ)il)f{x). (10) 



V 



/=1 



The index / labels different quark flavours, and the covariant derivative contains the gauge coupling g, 

D^ = {d^-igA^), A^ = T"A^(x), a=l,...iV2-l, F^,{x) = '-[D,,D,] . (11) 

The thermodynamic limit is obtained by sending the spatial three- volume V ^ oo. The difference to 
the Euclidean path integral at T = describing vacuum physics is the compactified temporal direction 
with a radius defining the inverse temperature, 1/T. To ensure Bose/Einstein statistics for bosons and 
the Pauli principle for fermions, the path integral is to be evaluated with periodic and anti-periodic 
boundary conditions in the temporal direction for bosons and fermions, respectively. 



A^(r,x) = A^(r + ^,x), ^(r,x) = -^(r + i,x) . (12) 



The path integral for vacuum QCD in infinite four-volume is smoothly recovered from this expression 
for T ^ 0. 

The partition function depends on the external macroscopic parameters T,V,fif, as well as on the 
microscopic parameters like quark masses and the coupling constant. The conserved quark numbers 
corresponding to the chemical potentials fif are 

Qf = i'noipf ■ (13) 

We will consider mostly two and three flavours of quarks, and always take rriu = rrid. The case = rriu^d 
is then denoted by Nf = 3, while Nf = 2 + 1 implies ^ mu^d- When the heavier charm quark is 
taken into account as well, we have Nf = 2 + 1 + 1. 

2.2 Centre symmetry 

Consider the partition function of pure gauge theory. The action is invariant under gauge transforma- 
tions g{x) which obey the same periodic boundary conditions, gij^y.) = gij + 1/T, x), 

K^^) = 9{^) [a, + '-d,^ gKx) ■ (14) 

It is also invariant under gauge transformations g' , which are periodic up to a global SU {N) matrix h, 
g'{T,x.) = hg{T + 1/T, x). The gauge transformed field satisfies the boundary condition, 

A^'(r + l/T,x) = /iA^'(r,x)/it, (15) 
5 



which is the originally demanded one if G Z{N) is in the centre of the group, h = zl = exp{i27rk/N), 
k e {0, l...A^ — 1}. Hence pure gauge theory allows for gauge transformations which are periodic up to 
a twist factor from the centre of the group. For fermions one has instead 

(r + 1/T, x) = -zips' x) , (16) 

which is anti-periodic for z = 1 only, i.e. fermions break the centre symmetry of the pure gauge action. 
The order parameter for centre symmetry is the Polyakov loop, i.e. a traced temporal Wilson line, 

L(x) = Triy(x) = T'exp (^g ' dTAi{x)^ , (17) 

which corresponds to the propagator of a static quark in euclidean time wrapping around the torus. 
Under gauge transformations, W^^(x) — g'(x)Vl^(x)g'~^(x) it picks up a centre element when transformed 
with a winding transformation, 

= L, L^' = z*L . (18) 

The expectation value of the Polyakov loop gives the free energy difference between a Yang-Mills plasma 
with and without the static quark. 

(L) ^- [ DU TrW e-^o = ^ = e-(^«-^°)/^ , (19) 
Z J Z 

For T = pure gauge theory is confining and it costs infinite energy to remove the quark to infinity, 
i.e. Fq = oo and therefore (L) = 0. For T ^ oo one has (L) ^ 0, which is no longer invariant under 
centre transformations and signals the spontaneous breaking of centre symmetry. Therefore, QCD in the 
quenched limit has a true (non-analytic) deconfinement phase transition corresponding to the breaking 
of the global centre symmetry. When finite mass fermions are added, the centre symmetry is broken 
explicitly, i.e. (L) ^ always and there is no need for a non-analytic phase transition. 

2.3 Chiral symmetry 

The QCD action is invariant under U{1) transformations on the quark fields 

■0^ — exp(iQ;) i/jf , (20) 

the corresponding conserved charge is baryon number. For n/ degenerate quark fiavours, the corre- 
sponding spinors can be written as an n/-plet ip = {ipf,...ipnf), and the QCD action is symmetric 
under the global SU{nf) vector rotations 

ip'f, ^ex.p{ie''T'')f,f ijjf , a = l,...n5-l, . (21) 

In nature, ~ and this is SU{2) isospin symmetry. For Uf zero mass flavours QCD is furthermore 
invariant under the corresponding C/a(1) and SUA{nf) axial transformations, 

^} = exp(m75) iPf , = exp(^rT»75) , a = 1, . . . nJ - 1 . (22) 

Hence, in the chiral limit the total symmetry is 

C/(1)b X Uil)A X SUinf)L X SU^r , (23) 

where the L, R denote the left and right handed linear combinations of the vector and axial-vector 
rotations. The axial U{1)a is anomalous, quantum corrections break it down to Z{nf). The remainder 
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Figure 1: Left: Schematic phase transition behaviour of Nf = 2 + 1 QCD for different choices of quark 
masses at fi = 0. On finer lattices, the chiral critical line moves towards smaller quark masses. Right: 
Chiral and deconfinement critical surfaces swept by the critical lines as fi is turned on. Depending on 
its curvature, a QCD chiral critical point is present or absent. 

gets spontaneously broken to the diagonal subgroup, SU^Nf)^ x SU{Nf)ji — > SU{Nf)Y, giving rise to 
nj — I massless Goldstone bosons, the pions. 

At finite temperatures, chiral symmetry gets restored above some critical temperature Tc. The order 
parameter signalling chiral symmetry is the chiral condensate. 

It is nonzero for T < T^, when chiral symmetry is spontaneously broken, and zero for T > Tc. Hence, 
in the chiral limit there is a non-analytic finite temperature phase transition corresponding to chiral 
symmetry restoration. This is of first order for Nj > 3 and either second order with 0(4) universality 
or first order for Nf = 2, depending on the strength of the U{1)a anomaly at For non-zero quark 

masses, chiral symmetry is broken explicitly and the chiral condensate (ipipf) 7^ for all temperatures. 
Again, in this case there is no need for a non-analytic phase transition. 

2.4 The finite temperature phase transition of QCD 

The order of the finite temperature phase transition at zero density depends on the quark masses 
and is schematically shown in Fig. [1] (left). In the limits of zero and infinite quark masses (lower left 
and upper right corners), order parameters corresponding to the breaking of a global symmetry exist, 
and for three degenerate quarks one numerically finds first order phase transitions at small and large 
quark masses at some finite temperatures Tc{m). On the other hand, one observes an analytic crossover 
(without singularities in the thermodynamic functions) at intermediate quark masses, with second order 
boundary lines separating these regions. Both lines have been shown to belong to the Z{2) universality 
class of the 3d Ising model E]. The critical lines bound the quark mass regions featuring a chiral or 
deconfinement phase transition, and are called chiral and deconfinement critical lines, respectively. 

QCD with physical values for the quark masses breaks both the chiral and the centre symmetries 
explicitly. Hence, the finite temperature transition does not break or restore either symmetry and there 
is no reason for it to be non-analytical. Indeed, physical QCD is in the crossover region as established 
using the staggered fermion discretisation in two independent ways: i) by showing that the finite 
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temperature transition at the physical point is non-critical and analytic in the continuum limit [7], ii) 
by explicitly mapping out the chiral critical line (so far on coarse lattices only) ^8j. At finite density the 
situation is much less clear, since straightforward Monte Carlo simulations are impossible due to the 
sign problem, cf. Sec. [TOl For larger chemical potentials some techniques signal the presence of a critical 
point. However, these are not yet independently confirmed. On the other hand, for small /i/T< 1 all 
independent methods find that both the chiral and deconfinement transitions weaken with chemical 
potential, so there is no chiral critical point at small chemical potentials. Reviews of the current status 
can be found in P . 

Associated with the deconfinement and chiral phase transitions is a critical temperature Tc{mf,fif), 
which depends on the number of fiavours and quark masses. The method to locate a transition in 
statistical mechanics is to look for infiection points in the changes of observables 0(x), O G {L, iptpf, ■ ■ ■}, 
or for the peak in the associated susceptibilities. Note that for an analytic crossover the pseudo-critical 
couplings defined from different observables do not need to coincide. The partition function is analytic 
everywhere and there is no uniquely specified "transition" . This holds in particular for pseudo-critical 
couplings extracted from finite lattices. As the thermodynamic limit is approached, the couplings 
defined in different ways will merge where there is a non-analytic phase transition, and stay separate in 
the case of a crossover. 

2.5 Ideal gases: fermions, bosons and hadron resonances 

Much of our intuition for QCD at finite temperatures is based on the behaviour of free gases. For 
example, the qualitative expectation of a rise in the pressure as a signal for deconfinement is based on 
the ideal hadron resonance gas in the low temperature phase and an ideal gas of quarks and gluons in the 
deconfined phase. This is also true on the lattice. At high temperature perturbative investigations allow 
for a controlled analytic understanding of the leading cut-off effects on the equation of state. At low 
temperatures a strong coupling expansion provides new qualitative insights by showing the emergence 
of the hadron resonance gas as a true effective theory of QCD in the strong coupling limit. 
The partition function for an ideal gas reads 

lnZl{V,T)=r^u.y J ^ ln(l + r^e-^^-'^-)/^) , (25) 

where rj = —1 for bosons and 77 = 1 for fermions, Ui gives the number of degrees of freedom for the 
particle i and is the chemical potential for a conserved charge associated with the particle. In the 
massless limit the momentum integration can be done in closed form and gives the Stefan-Boltzmann 
pressure 

vr^ ( 1 bosons , . 

V = Vi — T x< 7 r ■ ■ 26 
^ 90 \ { fermions ^ ' 

For an ideal gas of gluons and massless quarks we have two polarisation states for each of the (A^^ — 1) 
gluons, two polarisation states per quark, colours per quark and the same for anti-quarks. 

This is the Stefan-Boltzmann limit of the QCD pressure which is valid for vanishing coupling (7 = 0, 
i.e. in the limit T — )■ 00. 

There is another ideal gas system that is useful to model the low temperature behaviour of QCD, 
assuming that at low temperatures, when we still have confinement, the QCD partition function is close 
to that of a free gas of hadrons. While strong coupling effects are responsible for confinement of the 
quarks and gluons, the colour forces are screened beyond the hadron radius and interactions between 
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colour neutral hadrons are considerably weaker. Neglecting them, the QCD partition function factorises 
into one-particle partition functions, In Z{V,T, fii) ^ In Z/(V, T, /ij), and the pressure reads 

P = ^ Yl ln^.'(^,T,/x,) + ^ Yl (28) 

iG{mesons} i£{baryons} 

with 1] = —1,1 in the mesonic and baryonic terms, respectively. The particles to be inserted are the 
known hadrons and hadron resonances (for QCD purposes electroweak decays are to be neglected). 
Taking experimental values for the energies from the particle data group [H], one can supply this 
formula with hundreds of hadron resonances, do the momentum integral and obtain a thermodynamic 
pressure that compares remarkably well with Monte Carlo simulations of QCD [12]. It was later shown 
that this is no accident, but that the hadron resonance gas model can actually be derived from lattice 
QCD as an effective theory for low temperatures, i.e. strong coupling ^I3], see Sec. I3.12[ 

3 Lattice QCD at finite temperature 

Experience has shown that in order to avoid misleading conclusions, it is most important to understand 
and control systematic errors due to the lattice discretisation when attempting to extract continuum 
physics from results of numerical simulations. For this reason it is worth discussing the discretisation 
in some detail. 

3.1 Lattice gauge action 

We consider the lattice formulation of a SU {N) pure gauge theory on a hypercubic lattice, A^^ x A^,-, with 
lattice spacing a and spatial and temporal lengths L = aNs,T~^ = aNr, respectively. The standard 
Wilson gauge action reads 

SAU] = J2 E /3 (l - ^ReTrt/^,(x)V (29) 

where U^y{x) = Up{x) = U^{x)Uy{x + afi)Uj^{x + aP)Ul{x) is the elementary plaquette, and the bare 
lattice and continuum gauge couplings are related by /3 = 2N/ g"^. For numerical simulations one imposes 
periodic boundary conditions in all directions, 

f/^(r, x) = U^{t + iV,, x), [/^(r, x) = U,,{r, x + iV,). (30) 

To see the relation to the continuum action, write the links as infinitesimal parallel transport between 
neighbouring lattice points. Inserting this in the plaquette and expanding down the exponential yields, 

U^{x) = exp{-iagAf,{x)) , f/^^(x) = exp (za^F^j,(a;) + 0{a^)) . (31) 

In the limit of small lattice spacings the action reproduces the classical continuum action, 

^^[f^] = ^ E E (Tr^i(^) + ^(«')) • (32) 

X fl,U 

However, there are infinitely many more terms suppressed by increasing powers of the lattice spacing. 
As long as this is finite, the action and all results differ from the continuum action by these discretisation 
effects. As a consequence, the lattice action is not unique, only its classical continuum limit is. One 
may add or subtract any terms that vanish in the continuum limit without affecting the continuum 
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physics. This is made use of in the construction of improved lattice actions, which attempt to reduce 
discretisation effects. 

The connection between zero and finite temperature physics is most easily exhibited by the transfer 
matrix, which evolves quantum mechanical states by one lattice spacing in Euclidean time. It relates 
the path integral representation of a Euclidean lattice field theory to the Hamiltonian formulation and 
is completely specified by the lattice action [H], to put the partition function in its quantum mechanical 
form 

Z = f r(T^-) = f r(e-^-"^), ^ = aN, . (33) 

In this form we can immediately see that Z is equivalent to the partition function of a thermal system 
if we identify the temporal lattice extent with inverse temperature. The thermal expectation value of 
an observable is then 

(O) = Z-tr(e-f O) = $:(n|T-^0|n) = ^"i^'^'!];!"^"" • (34) 

As in the continuum, we are interested in the thermodynamic limit and hence Ng ^ oo, while keeping 
qNt = 1/T finite. In this form we easily see the connection to T = physics: projection on the vacuum 
expectation value is achieved by taking A^,- to infinity, with Eq the vacuum energy, 

V (n\0\n) e~°'^^^^"~^'^'> 
{0\O\0)= hm 2.nW^I^;e _ ^35) 

In order to describe our gauge theory at finite temperatures, we simply need to dispense with this step. 
The spectral decomposition of a thermal expectation value in Eq. (1341) gives some important qualitative 
insights. Its value is governed by the eigenvalues and eigenstates of the Hamiltonian, i.e. the spectrum 
of the theory at zero temperature. The effect of temperature is the Boltzmann-weighted contribution 
of higher states. This makes it plausible why for T < Tc the hadron resonance gas should give a good 
description of thermal expectation values: the En in this case are the masses, momentum and scattering 
eigenstates of hadrons, and the hadron resonance gas description only neglects the latter. 

This argument is valid up to the phase transition. It can be extended by defining a transfer matrix 
and Hamiltonian in a spatial direction, Z = Tr(e~"^^^^). is now hidden in the definition of the 
Hamiltonian. Thermal physics is thus equivalent to the "zero temperature" [N^ — ?■ oo) physics of the 
Hamiltonian H^, which acts on states defined on a space with two infinite and one finite, compactified 
direction. For large N^- the eigenstates of are equivalent to those of H, i.e. the vacuum hadron states. 
Now finite temperature is nothing but a finite size effect: thermal effects become noticeable once Nr is 
small enough for the system to be sensitive to the boundary. For T > Tc, the eigenstates of Hz are the 
screening masses m T. These general considerations remain the same when fermions are included. 

3.2 Including fermions 

Once a suitable fermion action 

^S^^Y^ ^/(a;)^xs/(m/) ipfiy) (36) 
f 

has been selected and minding the appropriate anti-periodic boundary conditions in the temporal di- 
rection, the Gauss integral can be done and we end up with the partition function 



Z{Ns,Nr;(3,mf) = ^ Dt/ ]^ e"^«[^l~^^[^''^^''^-f] = ^ Dt/ JJ det M(m/) e" 



Sa[U] 



f " f 

t/^(r,x) = U^ir + Nr,x), V^(r, x) = -^^(r + A^,, x) . (37) 
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However, the selection of an appropriate lattice action is far more intricate than in the pure gauge 
case. Even more difficult is the treatment of finite baryon chemical potential. Here we set /x^ = and 
postpone the discussion of finite density to Sec. [101 



3.3 Naive lattice fermions 

In the naive discretisation, one simply replaces the partial derivative with a symmetrised finite difference, 

V,/(x) = ^^^ + ^^~/^^~^^ (d, + laX + Oia^)^fix). (38) 

which introduces leading order (9(a^)-corrections to the action. Making the derivative covariant, the 
Dirac operator reads (with U_^{x) = WJx — fi)) 



/ ' \ ' 1 

Dxy{m) = ^ 7MV^(f/) + m = ^ 7^— (f/^(a;)5:r+^,y - U-f,{x)6x-ix,y) + m 5xy 



(39) 



The corresponding Fourier transform and quark propagator are 

. 4 

D{p) = m + -^7^sin(ap/,), 

X ^ m-m-^^^7^sin(ap^) ^ m - ig-^ J2^,1^^P^^ 
m? + a~2 sin^(ap^) + a~2p2 

with the periodic lattice momentum = sin(ap^), which is a consequence of the discrete Fourier 
transform. Consequently, in the chiral limit m = 0, the fermion propagator has 16 poles {ap^ = {0, vr} 
per space-time direction /i). It thus describes 16 degenerate fiavours of massless quarks where only 
one is desired. Since this is true for any lattice spacing however small, those 16 fiavours survive an 
extrapolation to zero lattice spacing, thus leading to an incorrect continuum limit with an enlarged 
particle content. The same observation holds for m 7^ with poles shifted away from zero. There are 
several ways of dealing with this problem, see e.g. IT6] . In the following we shall focus on Wilson 
fermions and staggered fermions, since they allow for the most realistic numerical simulations presently. 
For all fermion species, at finite temperature the temporal component of momentum is represented by 
the Matsubara frequencies Un = (2n + l)7rT, refiecting the anti-periodic boundary conditions for the 
spinor fields in the temperature direction. 

3.4 Wilson fermions 

A straightforward cure of the doubling problem is due to Wilson. His lattice Dirac operator is given by 

±4 

DZim) =lm+^) 6.,y - V (1 - ^,)U,{x)6.^f,,y , (41) 



(4 \ 1 



with the notation 7_^ = — 7/^. Note the additional mass term ~ 4/a. The corresponding Fourier 
transform and quark propagator read 

D^{p) = m + - 7^ sin(ap^) + - V'(l - cos(ap^)) 



p2 + {\p^ + am)' 



[D^)-\p) = — . "^r. 'r \ 2 — , (42) 
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with Pn = 2 sin(ap^/2). The third term in the Fourier transform is the Wilson term, which vanishes for 
p = but contributes 2/ a whenever a momentum component = vr/a. Hence, the mass degeneracy 
of the doublers is hfted. In the massless hmit there is only one pole at = while the 15 others have 
poles ~ a~^, with appropriate shifts for m ^ 0. So for a — )■ the doublers become infinitely heavy 
and decouple from the theory. Therefore, on sufficiently fine lattices Wilson fermions have no doubling 
problem. 

The prize to pay is a complete loss of chiral symmetry, since the Wilson term is not invariant under 
the transformations Eq. (1221) . That is, even in the massless limit there is no chiral symmetry and 
hence no spontaneous breaking of it either. As a consequence, the quark mass receives both additive 
and multiplicative renormalisation. One still may obtain massless pions by an appropriate tuning of 
the bare parameters, but they do not correspond to Goldstone bosons for any finite lattice spacing. 
Consequently, there is no true chiral symmetry restoration at finite temperatures. Acting as a mass 
and hence symmetry breaking term, the Wilson term should soften any chiral phase transition, similar 
to the quark masses in Fig. [TJ Fortunately, since the Wilson term decouples in the continuum limit, all 
continuum symmetries get gradually restored as the lattice spacing vanishes. One would then expect 
that Wilson fermions are more cumbersome when it comes to physics related to chiral symmetry, but 
that eventually the correct continuum physics is obtained. 

Finally, for numerical simulations it is often more convenient to express the Wilson Dirac operator 
in a slightly different form by a rescaling of fields and introducing the hopping parameter k, 

DZ = il-KH,y), «: = ^— L-^, H,y = J25y,.+f.{l + l,)U,{x) . (43) 
3.5 Staggered fermions 

The staggered action is obtained from the naive action Eq. (|39l) by a series of field transformations 
mixing lattice and Dirac indices of the original spinor fields, for detailed formulae see [IHl [16]. The 
result is diagonal in the spinor indices, i.e. corresponds to four identical terms in the action for each 
value of the spinor index. The fermion doubling problem may then be reduced by simply discarding 
three of those terms. The resulting action can be written as 

Sf = J2x{x)D,yMx{y), (44) 

with the staggered Dirac operator 

D%im) = -'^ri^iS^+f.^yU^ix) - 5^,y+/if/^(y)^) + S^yam , (45) 
and the staggered sign function 

ri^ix) = ni{x) = 1 . (46) 

Note that this Dirac operator contains no 7-matrices and the fields x{^) the action have no Dirac 
indices, i.e. are one-component fields. Hence, this action is expected to describe only four fermion 
species compared to the 16 of the naive fermion action. 

As a by-product of this reduction, the one-component fermion field is significantly cheaper in nu- 
merical simulations than four-component spinors. A further advantage of the staggered discretisation is 
that it retains a part of the continuum chiral symmetry. Defining 775 (x) to be a staggered phase playing 
the role of the 75 matrix, 

ry5(x) = (_i)-i+-2+-3+-4 ^ (47) 
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one verifies that in the chiral hmit, m = 0, the action Eq. f l44j) is invariant under continuous transfor- 
mations with a G M, 

^(a;) e*"'?5W^(3,)^ ^(3,) ^ ^(a;)e-''5(3.) . (48) 

While Eq. (jUj) thus offers a few advantages, one must also make sure that its reduced degrees of 
freedom describe Dirac fermions as desired. In order to see this, let us simplify without loss of generality 
to the non- interacting case, i.e. U^{x) = 1. The action can be reformulated once more by transforming 

the fields back into Dirac spinors ipa^ p3| [16], with a spinor index a = 1, ... 4 and t = 1, ... 4 labelling 
different fermion species, 



S, = i,'^ 



Y, m >/'<'>(.-)>/'"'(2) + Yl ■iA"'(-)7,V^'i/.<'>(.-) 

.t=l \ fl=l 



b ' 



2 



t,t'=l At=l 



(49) 



These Dirac spinors live on a blocked lattice with twice the original lattice spacing, b = 2a, implying a 
new coordinate z^ = 0,l,...A^s/2 — 1, and finite differences over b, 



^ fi N f{z + fi)- f{z-fi) f{z + il)-2f{z) + f{z-fl) 
^^fiz) = ^ , AJ{z) = p , (50) 



where we have introduced matrices = 7^ , /i = 1, 2, . . . 5. The Fourier transform and quark propagator 
now read 

D{p) = m+ - ^7^sin(6p^) - -(1 - cos ( 6p^)) 75 r^rg 
^J.=l 

, ^ - ^b'^ J2l=i lf^sm{bpf,) + l-i^T^T^ siv?{bpj2) 

D [p) = J 2 . (51) 

fejE^=iSin (6p^/2) +m2 

Indeed, the poles at = vr/a have disappeared thanks to the factor 1/2 inside the sin-function of the 
denominator, and we are left with only one pole at = for each ■0^*^ The staggered action thus 
indeed describes A^^^ = 4 species of fermions. Since these are remnant doublers, they are usually called 
"tastes" in order to distinguish lattice artefacts from the physical fermion species called flavours. 

The kinetic and mass terms in the first line of Eq. (H9|) are diagonal and hence display a SU{4) taste 
symmetry. The third term breaks this symmetry by mixing different tastes, thus lifting the degeneracy 
of their masses. However, for a — )■ the last term is suppressed compared to the other two and vanishes. 
Hence, in the continuum limit the four tastes become mass-degenerate and identical. This motivates 
the final step: doing the Gaussian fermionic path integral, the staggered action gives a determinant 
representing four tastes. In order to describe one single flavour it is then customary to take the fourth 
root of the staggered determinant, whereas two mass-degenerate flavours are approximated by taking 
the square root. Hence, for Nf = 2 -|- 1 staggered fermions, 

J DijDij e"^/ [^'^'^] = det[D'\mud)]^/^ det[D^*(m,)]^/^ . (52) 

It must be stressed that this last step is only well- motivated in the continuum limit, when the tastes 
are degenerate. In a numerical simulation at finite lattice spacing, this is not the case and one has to 
worry about the systematic effect that this introduces. With four tastes of quarks, there are 16 meson 
states which fall into eight multiplets [T7], characterised by different masses and degeneracies. These 
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index 




multiplicity rij 





75 


1 


1 


7o75 


1 


2 


7i75 


3 


3 


7i7i 


3 


4 


7i7o 


3 


5 


7i 


3 


6 


7o 


1 


7 


1 


1 



Table 1: Pion taste multiplets labeled by taste matrices have multiplicity rij within the multiplets. 

can be ordered according to the generator matrices that distinguish between the tastes, F^, which can 
be chosen to be products of 7-matrices [T71 [18], as in Tabled! Of the 16 pions listed in the table, only 
the one corresponding to = 75 is a Goldstone boson, i.e. its mass rriG vanishes for vanishing quark 
mass. Of course, this is consistent with the remnant pattern of chiral symmetry breaking for staggered 
fermions, where only one generator gets broken. The other pseudo-scalar mesons are non-Goldstone 
bosons and have masses 

mf = tuq + Smf, with 5mf ~ aa^ . (53) 

On coarse lattices with lattice spacing a > 0.25 fm some can be as heavy as several hundred MeV, as 
we shall see in detailed numerical analyses in Sec. 14.91 

Aside from the taste splitting, there are several other conceptual questions concerning staggered 
ferminons which are still under investigation and debate. It is not entirely clear whether the rooted 
determinants can be represented by a local lattice field theory when re-exponentiated. This has reper- 
cussions on whether the continuum limit of the staggered theory falls into the correct universality class. 
Another question concerns the coupling between topological sectors and zero modes of the staggered 
Dirac operator, since the coordinate and spinor degrees of freedom are staggered on different lattice 
sites. This is important in the small mass region close to the chiral limit, when the U{1)a anomaly plays 
a role. For references discussing these issues, see [19l EOl EH |22]. Since so far there are no conflicting 
results between staggered and Wilson fermions, the working assumption for the purposes of this review 
is that these latter questions can be satisfactorily answered one day, leaving cut-off effects and taste 
breaking as systematic errors to monitor in numerical simulations. 

3.6 Observables related to centre and chiral symmetry on the lattice 

In Sec. 12.21 we have seen that the centre symmetry of the partition function corresponds to a subgroup 
of the gauge symmetry of Yang-Mills theory. Since gauge symmetry is fully realised in the lattice 
formulation, so is centre symmetry. The Polyakov loop is represented by a product of temporal links, 

L{^) = l[Uoix) , (54) 

XQ 

with the same transformation behaviour as in the continuum. Also in analogy to the continuum, centre 
symmetry on the lattice is broken by including finite mass fermions. While the Polyakov loop is no 
longer an order parameter for confinement in this situation, it nevertheless signals a change in the 
dynamics by measuring the free energy of a static quark in the plasma. 

The self energy of the static quarks cause a UV-divergence, and for comparisons between different 
actions one needs to consider renormalised quantities. Since the Polyakov loop is not directly related 
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to an experimentally observable quantity, there is some freedom in its renormalisation, 

Ln{T) = Z{(3f^m . (55) 

One possible choice to fix the renormalisation factor Z[(]) is to relate it to the free energy of a static quark 
anti-quark pair at infinite separation, Lr{T) = exp (— Foo(T)/2T) [23]. A more practical procedure is 
to normalise the zero temperature static potential, obtained from Wilson loops, at a given distance to 
a prescribed value, e.g. Vr^vo) = 0, and then use the same additive shift for the free energy obtained 
from Polyakov loops at the same lattice spacing [23]. This translates into Z{P) = exp(\^(ro, /3)/2), with 
the unrenormalised zero temperature static potential. 

Similarly, the continuum chiral symmetry is explicitly broken by the quark masses and the finite 
temperature transition is a smooth crossover instead of a phase transition. While the chiral condensate 
for a light fiavour is no longer a true order parameter, it still signals changes in the dynamics. It is 
usually expressed by the inverse propagator, 

1 1 

for the light and strange flavours, respectively. Other interesting fermionic quantities are the corre- 
sponding quark number susceptibilities, 

Xf _ 1 d'^lnZ 



(^^/) = 77F;n^Tr(D7^), / = /,., (56) 



T2 VT^ difif/Tf 



, / = /, s . (57) 



As we have already discussed, the lattice situation for chiral symmetry differs significantly from the 
continuum. For Wilson fermions, chiral symmetry is broken explicitly for any finite lattice spacing, even 
for massless quarks, whereas for Nf = 2 staggered fermions it is reduced to a 0(2) subgroup. Due to the 
breaking of chiral symmetry on the lattice, both additive and multiplicative renormalisation of the chiral 
condensate are necessary. Renormalisation of the chiral condensate for Wilson fermions is discussed in 
[25j . the presentation here follows [261127]. The dimension three operator ipil} features a divergence ~ 
which can be removed by subtracting the corresponding vacuum value. The multiplicative divergence 
gets removed by multiplying appropriate factors of quark masses displaying the same divergence, ensured 
by axial Ward identities. Near the continuum limit the following relations hold between renormalised 
and PC AC quark masses, 

mR{i}^)R{T) = mpcAcZA^ip^iT) + 0{a) (58) 
mR{^P^P)R{T) = 2NfmlcAcZl^PpiT) + 0ia) (59) 

where Z^ is a finite renormalisation constant. The latter can be eliminated together with mpcAc to 
arrive at 

For the condensate one has 

A^^(T) = (^^)(T)-(^^)(T = 0), 

App(r) = I d'x{Po{x)PoiO)){T) - [ d'x{Po{x)Pom{T = 0) , (61) 



where (ipip) and Po{x) are the bare chiral condensate and bare pseudo-scalar density, respectively [25j . 
Note that Eqs. fl59[ I^U]) have different corrections and hence different scaling behaviour, which must be 
determined numerically at this point. Finally, a definition for a finite condensate which is used in many 
analyses employing staggered fermions is [28] 

(Wi,o-^W).,o 
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3.7 Tuning temperature, continuum limit and lines of constant physics 

According to Eq. f l5^ . the simplest way of tuning temperature in a numerical simulation on the lattice 
is by choosing for a given lattice spacing a. This is known as the fixed scale approach whose 
application to the equation of state is discussed in Sec. I5.3[ Often this is not satisfactory as this allows 
only discrete temperatures. Then it is more practical to keep Nr fixed and vary the lattice spacing a 
via the coupling, /3 = 2N/g'^{a), thus affecting temperature. This implies a few marked differences to 
simulations at zero temperature. In particular, simulation points at different temperature correspond 
to different lattice spacings and thus have different cut-off effects. The continuum limit, a — )• 0, via the 
running coupling corresponds to the limit /3 — oo while keeping observables fixed in physical units. For 
simulations at fixed temperature T, this also implies A^^ — > oo. A continuum extrapolation therefore 
requires simulations on a sequence of lattices with different A^^- 

In order to take the continuum limit, we need to compute expectation values of our observables of 
interest, (0)(/3, mf), in the thermodynamic limit for various lattice spacings a. Then we can extrapolate 
to a — )■ while keeping temperature and physical masses fixed. Similarly, when simulating different 
temperatures on a lattice with fixed Nr the lattice spacing is different for different temperatures. This 
means that also the bare quark masses in lattice units, arrif, need to be tuned such that the masses 
in physical units stay constant. A line of constant physics is an implicitly defined curve in parameter 
space along which the bare quark masses are tuned together with the lattice gauge coupling, such that 
the renormalised quark masses in physical units is kept constant, 

mf{amu,d{l3)^ ams{l3)^ l3) = const . (63) 

Renormalised quark masses are scheme dependent and complicated to compute. In practice it is better 
to define the lines of constant physics by physical observables that can be straightforwardly computed 
in lattice simulations, such as hadron masses or decay constants. These are renormalisation group 
invariants and the corresponding lines of constant physics read 

Of'^^\amu4{f3),ams{f3), (3) = const, i = l,2,3, (64) 

where we need as many observables as parameters to tune. In principle, all chosen observables must lead 
to the same result. However, in practice the presence of cut-off effects also in those quantities renders 
different choices inequivalent. For finite lattices spacings, observables display lattice corrections, 

OPhy«(a) = OP'^y^ + cia + csa' + . . . , (65) 

with different coefficients for different observables (and actions). Keeping some observables fixed in 
physical units will have other observables move while changing the lattice spacing. Thus different lines 
of constant physics become equivalent only in the continuum limit, and careful choices have to be made 
in order to minimise cut-off effects during the extrapolation. 

The relation between the lattice spacing a and the bare parameters P,mf is given by the renormali- 
sation group. E.g., in units of the Lambda-parameter in lattice regularisation we have to leading order 
in perturbation theory for SU{3) pure gauge theory, 

-bi/2bl 



aAL = ( ^ 1 e 1260, 

^0 = -^(n-'^Nf], h ^ ^ 



167r2 V 3 ' / V 167r 



102 - f 10 + ^ ) A^/ 



(66) 



Similar relations exist for the running of the quark masses. However, perturbation theory is not 
convergent for most accessible lattice spacings. The way out is to express the calculated observables in 
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terms of known physical quantities of the same mass dimension. For example, if we want to compute 
the pseudo-critical temperature of the QCD phase transition, T^, by keeping A^,- fixed and tuning /3, the 
location of a phase transition will be given as a critical coupling Pc- The critical temperature in units 
of a hadron mass is then 

T 1 1 

-^c _ (Q7) 

tuh acTTinNr a{l3c)mHNr ' 

In order to set the physical scale for Tc, we then have to calculate the zero temperature hadron mass 
in lattice units at the value of the critical coupling, {amH){(ic)- Equating this number with the hadron 
mass in physical units times the lattice spacing provides the scale for the lattice, 

mj,[MeV] 

This is easily translated back to the critical temperature in physical units. Of course, the same procedure 
applies for observables other than the temperature, and besides a choice of hadron masses, also decay 
constants or other observables that are easy to measure on the lattice may be used. 

This procedure is good as long as simulations take place at physical quark masses. For many 
problems, this is not yet the case, and furthermore there are many interesting theoretical questions 
concerning regimes with unphysical quark masses, such as the quenched and the chiral limits. In order 
to set a scale in those cases, one uses quantities that display only little sensitivity to the quark mass 
values. Frequently used is the Sommer scale tq and its variant ri, which are based on the potential of 
a static quark anti-quark pair evaluated to prescribed values 



r^'!^^) =1.65, (r-^^^) =1.0. (69) 

/ r=ro V dr J ^^^^ 

The advantage of these quantities is that they can be straightforwardly computed in lattice calculations 
using Wilson loops, and that their numerical values change indeed very little with variations of the sea 
quark masses over wide ranges. The disadvantage is that the signal to noise ratio for the potential is 
worse than for typical hadron masses, and some care has to be taken concerning rotational symmetry 
when using improved estimators. Moreover, numerical determinations yield dimensionless numbers 
in lattice units, ro/a,ri/a, but there are no direct experimental measurements for these quantities. 
Therefore, additional lattice simulations are necessary in order to establish a relation to experimental 
observables, introducing additional systematics. The closest relation of the static potential is to spectral 
quantities of heavy quarkonia. Examples used in recent precision calculations [30] are the mass splitting 
of different radial excitations of bottomonium, the mass splitting between the Dg and rj^ mesons or the 
decay constant of the fictitious r^j^, which is in turn related to fn, fx- Another possibility is to set 
the scale for the lattice spacing by a simulation of a light hadron mass at physical quark mass values 
as described above, and then compute ro/a,ri/a in order to set the scale also away from the physical 
point. 

It should be clear that determining the lines of constant physics as well as setting the physical scale, 
albeit being performed at zero temperature, are an important part of a finite temperature simulation 
and play a major role in controlling and eliminating the cut-off dependence of numerical results. 



3.8 Constraints on lattice simulations 

Lattice simulations are beset by systematic errors and uncertainties, the two obvious being finite lattice 
spacing and finite volume. An important question before running a simulation or when interpreting its 
results is then: how large and fine does a lattice need to be for a particular problem? The Compton 
wave length of a hadron is proportional to its inverse mass , and the largest Compton wave length 
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constitutes the correlation length of the statistical system. To keep both finite size and discretisation 
errors reasonably small, we need to require 

a < mjj^ < aNs. (70) 

For the low T confined phase we thus need a lattice size of several inverse pion masses. The push to 
do physical quark masses is only just beginning to be a possibility on the most powerful machines and 
with the cheapest actions (i.e. staggered, with Wilson rapidly catching up). On the other hand, at high 
T screening masses scale as mn ~ T, thus 

AT-i < 1 < NsN;\ (71) 

While we desire large in order to have fine lattices for a given temperature, the spatial lattice 
size should be significantly larger than the temporal one. This combination of requirements limits the 
directly accessible temperatures to a few times Tc. 

To get some numbers, for a temperature T = 200 MeV ~ 1 fm~^, Eq. fl55]) with A^,- = 4, 8, 12 implies 
a lattice spacing a ~ 0.25, 0.125, 0.083 fm. We shall see later that the scaling region for the equation of 
state for most actions is entered only for a < 0.1 fm. Demanding only m^L > 2 means for the physical 
pion L > 2 • 1.5 fm. From this the required lattice sizes A^^ follows for A^^ = 10 to be at least Ng > 30. On 
the other hand, renormalisation requires also zero temperature calculations, i.e. as large as possible. 
For a = 0.1 fm, A''^ = 30 corresponds to T ~ 7 MeV. 



3.9 The ideal boson gas on the lattice 

Similar to the continuum, we can evaluate the ideal gas for a bosonic field on the lattice. Starting point 
is the equation employing the free propagator, 

InZo = --IndetA = -TrlnA-^ = y V j —^Ixiif + {amf) , (72) 

n=-Nr/2 a ^ ^ 

where now the integration is over the Brioullin zone and we have the lattice momenta, 

p2^4sin2(^)+4^sin2(^), Au^ = A J^^in^^) + (amy . (73) 

Due to the finite lattice spacing the Matsubara sum is only from —Nr/2, ....Nr/2 — 1. After performing 
the Matsubara sum (for details, see [31]), and a change of variables by cI; = sinh(ai?/2) the partition 
function reads 

In Zo = -V l [ ^ ln(l - e-^-"^) . (74) 

a 

Note that, despite the fact that the details of the calculation are quite different, this formula could have 
been obtained from the contiuum case, Eq. ( 125|) . by simply replacing momenta, energies and integration 
ranges by their lattice analogues. 

One can now study the approach to the continuum by expanding in small lattice spacing. The 
integration range of the coefficients then gets extended to infinity and the lattice corrections to the 
partition function are due to the lattice corrections to the dispersion relation E[p). The latter is 
defined by the zeroes of the denominator of the propagator. Writing 



E{p) = (p) + aE(i) + a^E^^) (p) + . . . , (p) = ^p' + m^, (75) 
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and expanding the dispersion relation 



sinh2(— ) = ^sin2(^)+^ 



2 



4 ' 

4+48 - 2^1, 4 48 4 +^^^^' ^^^^ 



one finds the leading a^-correction to the dispersion relation 



Note that this correction term breaks the continuum rotational invariance of the dispersion relation at 
the pf-level of momentum components. Changing to dimensionless variables, Xi = Pi/T,e = E/T, we 
have for the expansion of the pressure 



2 f d?x e^'^\x) 
T4 KT^Jcont J (27r)3e-(°'W - 1 + ' ' ' ^ 



P ( P 



Hence the pressure of a free gas of bosons has leading lattice corrections of 0{a?). For the massless case 
the integral can be done in closed form and one arrives at (for details, see [31] 



p 87r2 1 ^ / 1 , 

l + lT^l^ + OhrTi , (79) 



Pcont 21 iV2 

where Pcont is the continuum result Eq. ( |26|) . A similar calculation including finite size effects due to 
finite volume can be found in [Ml. 



3.10 The ideal fermion gas on the lattice 

This procedure can be generalised to the fermionic case [35j to obtain 



P_ Ar3 /"^ d^P 

2^4 



a 

For Wilson fermions, the inverse propagator is of the form 

D{pa,p) = + + arn^ , = sin(ap^) (81) 
Introducing again the energy variable sin^(ap4) = —sinh^ (ai?), the dispersion relation is 

(aE)2 + l(a^)4_^ f^^p^y _ (^\^amY-^ + -{aEYa''p'' + am {{aEf - a^p^) = , (82) 
12 Y 3 / 4 2 

i ^ ' 

which implies for the leading lattice correction 

= ^J^^\ ^^'^ = -^- (83) 
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Figure 2: Stefan-Boltzmann limit for massless fermions, normalized to the continuum result. The solid 
line is the analytic 0{1/N^) prediction. From [36] . 



The dispersion relation for staggered quarks to order pf is 

1 ^ 1 

{aE^ - a'p^ - {amy + - Y,im)^ + -(«^)' = (84) 

with the solution 

= -I {T.p' + E^'^'^ ■ (85) 

As in the pure gauge case, the corrections to the continuum result start at O(a^), whereas for Wilson 
fermions this is only true in the massless limit. The approach to the continuum between free massless 
staggered and Wilson fermions is shown in Fig. [2l Note that the scaling region, where leading lattice 
corrections are ~ = T'^/N^, sets in for A^^ > 10 only, for coarser lattices corrections seem to be smaller 
in absolute size for staggered fermions. 



3.11 The equation of state from strong coupling expansions 

Strong coupling expansions are well known from spin models and QCD at zero temperature. In contrast 
to the asymptotic series obtained by weak coupling expansions, they yield convergent series in the lattice 
gauge coupling /3 = 2A^/ g"^ within a finite radius of convergence. The series approximate the true answer 
the better the lower /3. For a fixed A^,- this means low temperatures and hence the confined phase. For 
pure gauge theory, of the deconfinement transition represents an upper bound on the radius of 
convergence. A detailed introduction to strong coupling methods in the vacuum can be found in [Ti] . 
Here we merely discuss modifications for finite temperature and the calculation of the pressure [371 US] ■ 
Using the formalism of moments and cumulants [H], the pressure of pure gauge theory on a lattice 
with temporal extent A^^ can be written as 

p(ArO = ^lnco(/3) + ^ J] a{C)l[<^{X,r, (86) 

c=(x"») « 

^In the first calculations [551 [33] the leading coefficient is incorrectly reported as 107r^/21. 
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Figure 3: Leading order tube contributing to the pressure for A^^ = 4. 



where C specifies a cluster of graphs and a{C) is a combinatorial factor [T^. is the contribution 

of a graph Xj, 

$(X,) = JdUH drCrXriUp) , (87) 



p€Xi 



with Xr{Up) the character of the d^-dimensional representation r of the plaquette Up, and Cr(/9) are 
expansion coefficients. The ones for the lowest- dimensional representations go as 



u 



/?2 

432 ^ 



and the fundamental representation coefficient Cf = u serves as an effective expansion parameter. The 
limit Nr ^ oo corrresponds to the zero temperature contribution. 
The physical pressure is given by the difference 



^9) 



where the subtraction of vacuum contributions is necessary for renormalisation. The lowest order graph 
Xi contributing to the difference is a tube of length Nr with a cross-section of one single plaquette, 
as shown in Fig. [3l There are three spatial directions for the cross section of the tube, giving a factor 
of 3. Translations in time take the graph into itself and do not give a new contribution, while we get 
V^{Xi) from all spatial translations, and a factor of two for the fundamental and anti- fundamental 
representation of SU{3). Thus, the leading order result for the pressure is 

If we now consider fixed Nr, low /3 and u{(3) corresponds to low temperature. We thus recognise the 
exponential suppression of the pressure in the low temperature, strong coupling regime. 

Corrections consist of plaquettes in higher dimensional representations inside the torus, additional 
cubes on the tube and slits by inserting two fundamental plaquettes at the same location. One finds 
for all Nr > 5 up to 0{u^) in the correction p!3] (with c = 1 + 3m + 6f + 8w — 18m^) 

3 



a p{Nr, u) 



N 



1 + 12 NrU^ + 42 Nru'' - Nru' 
^ ^ 2048 

597663 7 / ^ 2 72206061 > . 

NrVj + 83 Nr^ + Nr 1 



+b 



Nr 



2048 

1 + 12 NrU'^ + 30 NrU^ 



40960 



17191 
256 



NrU^ 



3819 

-180 XV + ( 83 Nr^ + Nr 1 



(91) 



Similar techniques can be used to include dynamical fermions p3]. Wilson's fermion action in the 
form of Eq. permits an expansion in powers of the hopping matrix, 

,J 



(92) 
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A truncation of this expansion is valid for heavy quarks, for which k/ is small. For sufficiently strong 
coupling, the gluonic contribution to the partition function can be omitted and 



Z{Kf) = / [dU] exp 



oo ; 

I f 1=1 



(93) 



In principle the sum extends over all closed loops I on the lattice, which for finite temperature lattices 
can have non-trivial winding numbers, resulting in Polyakov loops. Expanding all terms up to 0{k^^'') 
and doing the group integrals one finds for the pressure for two flavours up and down [13] , 



P = ^|4(2/s:„)2^^+8(2/c„2/c,)^-+4(2«:,)2^^| 
+ ^ |4(2/s:„)3^^ + 6 [{2Kuy2^d\ 

+ 6[2Ku{2K,)'f^ + 4(2«:,)3^-| (e^'^'^^- + e'^^'^^^-) . (94) 

3.12 The hadron resonance gas from strong coupling lattice QCD 

It is now interesting to ask how the QCD pressure can be interpreted in the strong coupling regime. It 
can be shown from flrst principles that it corresponds to a hadron resonance gas, without recourse to 
experimental numbers [13]. Since the hopping expansion is valid for heavy quarks, it is convenient to 
use the non-relativistic approximations E' ~ m + p^/ (2m). Inserting this in the expression for the ideal 
gas, the momentum integration gets trivial and the ideal gas pressure takes the form 

P. = jp^V In [1 + r^e-^"^'-^')^^] . (95) 
For non-relativistic particles it is also justifled to consider the limit m ^ T = l/N^-, which yields 

p = 5^p. = ^5^..e-(---)^^ (96) 

i i 

Using the fact that to each baryon with there is an anti-baryon with — /i^, we obtain 

P = ]^ E ^Me--^^ + ]^ E ^Be--^^ cosh(^^iV.), (97) 

M B 

where we have split the sum into a mesonic and a baryonic part. It is this form of the pressure which 
can be derived from a flrst principles strong coupling calculation. 

In the case of SU{?>) pure gauge theory, the three lowest lying glueball states, which can be extracted 
from the strong coupling series of plaquette correlators, have masses of [Ml |39] 

ami {A) = -4 In M - In (1 + 3m + 6vi + 8w - 18u^) + 0{u^) 
am2{E) = -4 In M - In (1 + 3m + Qvi + 8w - 18u^) + 0{u^) 
amsiT) = -41nM-ln(l-3M-6t;i + 8w) + C(M^). (98) 

The arguments of the masses denote the spin representation of the associated point group for lattice 
rotations, and the subscript the corresponding spin degeneracy of the glueball. The strong coupling 
series for the pressure, Eq. ( 19T|) . can now be rewritten as 

p = [e— + 2e-'"^(^)^^ + 3e-'"^(^)^^l (l + O (u^)) , (99) 
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which has to be compared with Eq. fl97j) . We see that the hadron resonance gas result equals the one 
derived with strong coupling methods to the orders considered here. Corrections to the ideal hadron 
gas results appear in O(m^). 

These considerations can be extended to include dynamical fermions. Using aruf ^ — ln2Kj, which 
is valid for heavy fermions |10lllT], one may derive the hadron masses to leading order hopping expansion 

Mesons: ^"^fp ~ — ln2Kj — ln2K// 

Baryons: amfffn = — ln2fi;j — ln2Kj/ — ln2Kj" . (100) 
The pressure, Eq. (194|) . may then be rewritten as 



P 



I 0- 1- J 

^i^Yl e-'"^^^)'^^ + 8 J2 e-"^^^^)"^^ I cosh (/x^iV.) , 

'^^ I 1+ 3+ J 



(101) 



where the sum over the pseudoscalar mesons includes the pions as well as the rj^, and the vector mesons 
include also the u^, since the calculation is valid for heavy quarks. The prefactors for the baryons 
include the spin degeneracy as well as a factor 2 counting the corresponding anti-baryons. 

We may thus conclude this section with the remarkable observation that, in the fully interacting 
theory, the hadron resonance gas emerges from a first principles calculation in the strong coupling limit of 
the equation of state for lattice QCD. Interactions between hadrons are suppressed by high powers of the 
inverse coupling l/g"^ until very close to the critical temperature. Conversely, agreement of simulation 
data with hadron resonance gas predictions signals that the underlying dynamics is governed by strong 
coupling. Furthermore, these results show from first principles that the pressure is exponentially small 
for temperatures not far below the quark hadron transition. 



4 Improved actions 
4.1 Generalities 

Since lattice actions are not unique, one may use the freedom in the choice of irrelevant terms to minimise 
cut-off effects. This can be systematised by making all a-dependence explicit. We are interested in a 
lattice action S^. Following Symanzik |12], one may write down an effective continuum action with the 
same symmetries as 5*^, 

^eff = J d^x {Coix) + aCi{x) + a'C2{x) + ...}, (102) 

where the Ck are composite local operators with canonical dimension A + k and Cq corresponds to 
the continuum theory. S^s thus produces the same vertex functions as Sl up to a certain order in a. 
Similarly, one may define effective fields as a power series in the cut-off, (/>eff(a;) = (pQ^x) + a0i(x) -|- . . . 
The idea of improvement is to construct lattice representatives of the 0^ and add them to the lattice 
action Sl and the fields in the observables of interest. If the coefficients are appropriately chosen, 
lattice corrections to a certain power in a can be made to disappear. The task of improvement lies in 
the choice of the coefficients. A straightforward way to determine them is lattice perturbation theory, 
but in this case the improvement is not complete as it holds to only a finite order in the coupling 
constant. Alternatively, one can determine the improvement coefficients non-perturbatively through 
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simulations upon which the improvement is exact to a given order in the lattice spacing. However, this 
is often complicated or associated with high computational cost. The only non-perturbatively complete 
improvements to date are (9(a)-improved Wilson fermion actions, see below. 

Besides the Symanzik program, there is also a more "empirical" way to improve lattice actions and 
operators. These are based on the idea that the discretisation effects arise through short wave length 
fluctuations of the gauge fields on the scale of the cut-off. Smoothing these fluctuations should decrease 
cut-off effects. This can be achieved by replacing elementary links by more extended objects, such as 
linear combinations of a link with the adjacent staples. These "smeared" or "fattened" objects have 
their dominant fluctuations on the scale of several lattice spacings instead of one. In some cases fat links 
are related to a perturbative improvement to some order in the coupling constant or lattice spacing, but 
in most others the improvement is not complete to any definite order in a. Rather, the improvement 
effect has to be judged experimentally by comparing calculations in different schemes. Besides the 
theoretical performance of a particular improvement scheme, its numerical efficiency is of course also 
an issue, but detailed comparative measurements in the thermodynamic context are missing to date. 
In the following sections, we discuss some schemes that have found applications in calculations of the 
equation of state. 



4.2 Improved gauge actions 

The Symanzik program for SU{N) gauge theories has been discussed in [331 Hi]- addition to the 
plaquette in the standard Wilson action Eq. (129|) . one may add larger planar loops of size x / in such a 
way that the leading 0{a?g^) corrections are subtracted and deviations from the continuum start with 
0{a^g^ ,a?g^). The improved gauge action then reads 



X l<fj,<l/ 



Co ( 1 - ^ReTrf/;f (x)] + c, (l - ^ReTrU'^:^' 



(103) 



where is the sum of the k x / and / x k rectangles. The coefficients in general are functions of the 
gauge coupling, Ci{g^), and in the weak coupling limit g'^ need to satisfy co(0) = 1 — A;^Z^Ci(0) in 
order to reproduce the correct classical continuum theory. For the choice 

the leading O(a^) corrections to the classical continuum limit are removed. Since this holds to leading 
order in perturbation theory, the action is said to be tree-level improved. The standard Symanzik 
improvement employs six-link rectangles, x / = 1 x 2. 

One can push this to higher orders in g'^ by adding further loops of larger length and non-planar 
loops. A non-perturbatively optimised expansion is obtained by multiplying the expansion coefficients 
with so called tadpole factors, which are to be determined self-consistently from plaquette expectation 
values [IS], 

ci,tad = u,'^'^^-'^cm, 4 = ( E - ^ReTrf/,.(x)') \ . (105) 

This amounts to a resummation of higher order terms to create a better expansion parameter, in the 
sense that convergence is faster. However, one should keep in mind that the resulting improvement is 
still not complete (i.e. not to all orders in the coupling constant) and depends on the efficiency of the 
resummation. Finally, renormalisation group methods have been used to design an action which is close 
to the true effective action in the renormalisation group sense along a renormalised line of constant 
physics [16]. This action can again be projected onto a two-parameter space as in Eq. f ll03p . 
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/3rt 


/5pg 


tree-level Symanzik 


5 
3 


1 

~6 





p4 


5 
3 


1 

~6 





asqtad 


1 


--^(1 + 0.4805a,) 


-0.1330^a, 


Iwasaki 


5 
3 


-0.662 






Table 2: Couplings x = pi, rt, pg defining the pure gauge part of the actions used in this review. 
Here Uq is the tadpole coefficient, Eq. fll05p . and a, = — 41n(uo)/3.0684. 



Often combinations of these possibilities are taken, and the gauge action can be written as [47] 



X,fl<l/ X,ll<l/<(7 



It consists of the ordinary plaquette as well as of planar and non-planar six-link terms, with fi, u, a E 
1, 2, 3, 4 and P^j^, i?^^ and C^i/o- denote the normalized trace of products of gauge field variables Ux^fj,- 



— J-Rp Tr ( TT TT - TT - - u'^ TT'^ TT^ 

-t-(/x ^ z/) + (z/ ^ a) + (/i 4-)- Z/, /i ^ -/i)) . (107) 

The choice of parameters for different versions used in the context of the equation of state are sum- 
marised in Table [2l 



4.3 Clover- improved Wilson fermions 

An improvement scheme for Wilson fermions which is frequently used in numerical simulations, both 
at zero and ffnite temperatures, are clover improved Wilson fermions [48] . 

Ssw = Sw + Csw iga'^ ^ ■^^{x)ai,uF^^ix)tp{x) , 

x,ii,u 

= ^ E (^M'^(^) - ' ^.^ = ^b.,7.]. (108) 

The additional term involves the clover of all four plaquettes sharing the lattice point x in each plane. 
In this case one uses the propagator of standard Wilson fermions, but with rescaled bare quark mass 
and gauge coupling [49] . 

= {m — mc){l + brnarn), g^-^ = g^il + hgam) . (109) 
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Here mc{g'^) denotes the additive quark mass shift corresponding to the chiral hmit, which hke the 
coefficients Cgw, bmig"^), bg{g^) can be determined order by order in perturbation theory, 

Csw = l + 0(^7'), b^ = -^ + Oig'), hg = 0{g^). (110) 
The tree-level dispersion relation at zero momentum reads 

E(0) = mil - ^am) + 0{a^) = m^"^ + 0{a^), (111) 

and is (9(a)-improved after rescaling the bare parameters. The tuning of Csw(/3) has also been done non- 
perturbatively, resulting in an interpolation formula which provides a completely (9(a)-improved 
Wilson fermion action in the /3-range covered. 



4.4 Twisted Mass Wilson fermions 

The twisted mass formulation is obtained by a chiral rotation in flavour space and requires an even 
number of pairwise degenerate flavours. For Nf = 2, after a rescaling of the fermion fields its action 
takes the form 

Stm = Sw + ifJ' i'l^T^i^- (112) 

Here /i denotes the twisted mass parameter and the diagonal Pauli matrix acts in flavour space. The 
bare quark mass ruq in this formulation is a combination of the bare mass m in the standard Wilson 
action and the twist parameter, 

rrig = a/tti^ + jj?, m = nig cos(co'), jj = nig sin(co') . (113) 

Twisted mass fermions provide an automatic (9(a)-improvement over Wilson fermions if the quark mass 
is solely determined by the twisted mass parameter /i, which is the case at maximal twist, u = tt/2. 
For a detailed derivation and recent review of twisted mass QCD, see [5T] . 
The free quark propagator for twisted mass fermions is 

Dtm{p) = — ^2 ' (114) 

P + + (^^) + (a/^)^ 

with the usual dimensionless lattice momenta 

= sin(ap^) , = 2sin(ap^/2) . (115) 

The standard Wilson propagator is recovered in the limit w — )■ 0. The dispersion relation E{p) can 
again be expanded in small lattice spacing and for vanishing spatial momentum we have 

E{0) = mg(l- ^amg cos(a;)) + 0{a^) . (116) 



2 

Cut-off effects of 0{a) set in with quark mass, and their removal at maximal twist, w = 7r/2, is apparent 



4.5 Naik and p4 improvement of staggered fermions 

Let us recall that for the staggered fermion action Eq. ( I44l) . the leading order corrections are (9(a^). 
Improvement of staggered fermions thus operates on the same level as that of the gauge action. Tree- 
level improvement for the staggered fermion action has been performed in [52] , by replacing the lattice 
covariant derivative by a higher order discretised version, cf. Eq. (l38l) . 

V^(f/) ^ V,iU) - ^Vj(f/) . (117) 
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- 5x-u+jv,y\ (118) 



The second term contains three hnks connecting fermion fields separated by one lattice spacing. A 
similar effect can be achieved when also including fermion fields separated by up to three lattice spacings. 
This can be illustrated denoting the free and massless staggered Dirac operator by 

H \i=l,3 v^^j=±2 

The constraints 

1 

Cl,0 + 3C3,o + 6Ci,2 = - 

ci,o + 3c3,o + 6ci,2 = 24ci,2 (119) 

ensure that the dispersion relation is rotationally invariant through order p^. For the Naik action, this 
is achieved by the choice ci^ = 9/16, ca^o = —1/48, with dispersion relation 



3 / . 1 



E = P+^\p'--X.P']^' + (^(^')^ (120) 







which is duly improved through 0{a ). The so-called p4 action corresponds to the choice Ci 
3/8,Ci^±2 = 1/48. While this particular choice still leaves a^-corrections in the action, it has the same 
dispersion relation Eq. fll20p . but smaller coefficients 0{a^). To improve the thermodynamics of ideal 
gases improved dispersion relations are sufficient, while higher order corretions are smaller than in the 
case of the Naik action. It should be stressed, however, that this construction relies on weak coupling 
which is realised only at very high temperatures. At lower temperature the perturbatively determined 
coefficients get increasingly out of tune with the renormalised couplings in the true effective action. 
Moreover, both the Naik and the p4 improvement have little effect on the taste splitting caused by the 
staggered action. 



4.6 Fat links for p4, asqtad and HISQ staggered actions 

Taste splitting at order can be cancelled by adding four-quark operators to the staggered action. 
However, a similar effect can be achieved by replacing the link in the lattice covariant derivative with 
specific fat links [M| The simplest fat links are obtained by replacing the link with itself plus the 
surrounding staples [56], such as the first term in the square brackets of Eq. f ll22p . Expanding down 
the gauge fields from the link variables, the fattened link adds a term [57j 

a^WstpD.F.^tP + 0{a^) (121) 

to the fermionic action in continuum notation (as well as to the staggered action after field transfor- 
mations). Clearly, it affects (9(a^)-improvement. From the pure gauge structure of this term, it will 
generally mix with similar terms appearing in improvement schemes for the pure gauge action, and care 
has to be taken as to appropriate combinations of the two. The staples used above are paths using 
three links. In a further extension also longer and non-planar paths are used [581 159] , 



'122) 



with the three-, five- and seven link terms. Fig. HJ 

Sj!J{x) = U,ix)U,ix + u)Ulix + ix) 
SUlix) = t/.(x)S(5(x + z>)f/t(x + A) 
S^JU^) = U^ix)S;^J^{x + u)Ul{x + ft) 
SitKx) = U,ix)S;^Jix + u)Ulix + fi). (123) 
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Figure 4: The simple link, three link staple, five link staple, seven link staple and double staple used 
in suppressing taste symmetry breaking. From [60] . 





Cl 


C3 


Cl2 


p4 


CO 1 00 





1 

96 


asqtad 


1 


1 





2 


48^2 



Table 3: Coefficients Cj for the fermionic p4 and asqtad actions, Eq. fll24p . as used in Sec. 17.31 



It can be shown that for appropriately tuned parameters with the first three terms the tree-level 
0{a'^) taste splitting is removed, but an a^-correction due to low momentum modes gets introduced by 
the seven link path. This is removed again by the double length staple S^^^ [61j. A staggered action 
including the Naik term, fat-5, fat-7 links and the double staple removes all tree-level errors, resulting 
in leading corretions of 0{a^^g'^a^) [GO]. The corresponding action is hence called asq, while the version 
with additionally tadpole improved coefficients is called asqtad. All these features can be summarised 
following the compact unified presentation (with changed normalisation of the coefficients) [17] 



m 5^y + (ci A[U]^y + C3 Bi[Ul^y + Ci2 B. 

^m(^) {yx,^lUx+^A.,^JLUx+2^A.,^l^x y~3fi - Ul_f^ JJl_2f,^^Ul_^^^^8x 



(124) 



11 UytfJ. 



Ux,^Ux+[i,uUx+ji+i>,u8x y-fi-2u Ux-u,u^x-2u,u^x-fi-2u,fj,^x y+fi+2v 



^x y+jl+2v 



+ {Ux,uUx+u,uUx+20,fj.8x y-fi-2u " f^i-/i,;,^i-/i_p,i.^i-/i-2i>,i.^: 
+ (^i-0,u^x-2^>,u^x~'2^',^J■^^ y-fi+2u ~ Ul_i^^^Ux^fi^uUx-jl+u,u8x y+fi-2u 
+ (^Ux,^Ul^~_jy^^Ul^~_20^^8x y-fl+20 — Ux,uUx+u,uUl_i^^2i>,fj.^x y+fi-2P 



A further significant improvement of the taste symmetry is achieved by the Highly Improved Stag- 
gered Quark action. In its original design in [62] it includes additional one-loop taste changing operators 
which are particularly useful for the spectroscopy of charmed mesons. For its use in thermodynamical 
simulations discussed here, it contains the same ingredients as the asqtad fermion action, but the fat 
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xw 



Figure 5: Left: Fat links and APE smearing replacing a link by itself and the sum of surrounding 
staples. From [57]. Right: Expansion to first order in the p^j^^ of the stout link variable U^^^ in terms of 
the original links. Closed loops include a trace. From 



links are reunitarised. (Note that the sum of link matrices is no longer unitary). Its use in conjunction 
with the tree-level Symanzik improved gauge action has been dubbed HISQ/tree action 



4.7 Stout links 

A more heuristic way to use smeared links in the Dirac operator is based on observations made in 
spectrum calculations at zero temperature, where smeared observables have less fluctuations on the 
scale on the cut-off. This suggests to use similar methods also on the level of lattice actions. One 
would expect that fluctuations on the scale of a lattice spacing get already damped in the generation 
of configurations, which is tantamount to reducing cut-off effects and hence improvement. 

The simplest possibility is called fuzzing or APE-smearing and consists of replacing a spatial link 
with the sum of itself and a weight factor times the sum of the adjacent staples, as in Fig. [5] (left) 
[56] . To maintain all symmetries of the original link, a projection back to SU{3) is applied. Empirical 
investigations show that the projection back to SU{3) is very desirable for the improvement of taste 
symmetry in staggered simulations. On the other hand, the projection step is non-analytic whereas for 
use in a Monte Carlo algorithm one would like to have different iable objects such that the response to 
infinitesimal changes can be computed. An analytic way to smear link variables goes under the name 
of stout links [64J. Consider the weighted sum of staples between neighbouring lattice sites. 




(125) 



The weights p^jy are tunable real parameters, with a common choice p = PikiP^iA = Pi/i = 0. Then an 
SU (3) matrix can be constructed in the following way 

fi^(x) = C^{x) Uj^{x), (no summation) , 

Q,{x) = '-(qI{x)-Q,{x))-^Tt(qI{x)~Q,{x)) . (126) 
This is the basis for an iterative smearing procedure, with a step from level n to ri -|- 1, 

Uj:^'\x) = exp(^Q(r)(x)) Ujr\x) . (127) 
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One can show that this step maintains all the relevant symmetries of the originial link. The first such 
smearing step is shown schematically in Fig. O (right). For a related smearing procedure constructing 
hypercubic links, see [05]. It should be stressed, though, that these methods do not implement a 
subtraction of terms of a specific power of the coupling or lattice spacing from the action, nor do we 
know a relation between smearing steps and the renormalisation group. Rather, their improvement 
effect has to be monitored in the results of simulations. 

4.8 Comparison of the weakly interacting gas for improved actions 

As the previous sections have illustrated, constructing and choosing improved actions is a complicated 
matter with many aspects to be taken into account. Let us now see how the different choices perform 
under different conditions. First, let us discuss the high temperature behaviour, which is affected by 
the perturbative improvements. The lattice pressure for free quarks can be expanded as a power series 
in the lattice spacing a, which at finite temperature is equivalent to expanding in l/N^, 



remarkably are the same for standard staggered and Wilson-type fermions. Differences are introduced 
only in higher orders of the lattice spacing. In particular, for the Naik and p4 improved staggered 
actions, corrections now start as 



A comparison of the free energy density of a gas of free massless staggered quarks is shown in Fig. [6] 
(left), where the standard, Naik and p4 actions are each normalised by the continuum result. Note 
that the leading cut-off effects of the Naik and p4 actions have the opposite sign compared to the 
standard action. It is apparent that the continuum value of the free energy density is approached more 
rapidly for the improved actions than for the bare action. However, one also sees that the improvement 
effect becomes significant only for sufficiently fine lattices, N^. > 6. Fig. [6] (right) shows the continuum 
approach between the p4 action and various actions using fat and smeared links. By design the p4, 
asqtad and HISQ actions remove O(a^) discretisation errors in the dispersion relation and hence in 
the thermodynamic functions, whereas stout smearing still shows a remnant a^-dependence in the free 
fermion gas. 

However, for a controlled continuum extrapolation, it is not necessarily important that cut-off effects 
are small in absolute size, but that their functional behaviour is controlled. The expected a^-behaviour 
is also shown by the stout action for > 10, allowing for a continuum extrapolation, whereas for 
coarser lattices a"^- correct ions become important. 

The ideal gas and massless limits are somewhat extreme simplifications and far from the cases of 
practical interest. In ^36j the effects of mass and interactions were taken into account. Note that in order 
to compare different discretisations in the massive case, the renormalised mass needs to be the same 
between different discretisation schemes. For non-vanishing masses the lattice corrections of Wilson 
actions to the free fermion gas start at 0{a), in contrast to the staggered actions. However, the linear 
coefficient is found to be exceedingly small compared to the O(a^) contribution which dominates in 
Fig. [7] (left). This also follows from the fact that the C(a)-improved Wilson actions do not scale better 




For massless quarks, the coefficients p^^-^^T^ and p^'^'^^T'^ can be calculated in closed form [35] and 




(129) 
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Figure 6: The free energy density of an ideal quark gas calculated for different values of the temporal 
extent Nt- divided by the corresponding value for iV^ = oo. From ^53j (left), [63j (right). 



than the unimproved one. Thus, for the ideal gas Wilson actions scale comparably to the unimproved 
staggered action and a continuum extrapolation using a^-behaviour is possible for A^,- > 10. 

The leading 0{g'^) corrections to the fermionic pressure due to interactions are given by two diagrams. 

The result for massive fermions is shown in Fig. [7] (right). Despite the rather large numerical uncer- 
tainties, one finds a qualitative correspondence to non-interacting case, especially the small dependence 
on the quark mass and the particular discretisation. As in the free case, on lattices with Nr > 8 
Wilson fermions are found to be competitive with unimproved staggered fermions. Again, scaling 
is not observed for any of the discretizations before > 10. Note that for Nr < 8 the absolute size 
of the cutoff effects for the unimproved Wilson action appears smaller than for the clover action. This 
indicates significant higher order contributions for the A'^^ considered here. 

The comparison of the scales between the free and the interacting situation shows that considerations 
based on the ideal gas limit alone may be insufficient for the discussion of improvement schemes. In a 
regime, where perturbation theory is valid, we can consider the relative cutoff error for the pressure, 

\^Pf \ /PF,cont = \PFM - PF,cont \ /PF.cont, (131) 

where PF,cont and p^^iat are each calculated to some specified order in the coupling. In the weak coupling 
limit the denominator is dominated by the ideal gas limit, and the cutoff effects of the interactions 
get normalised to that number. However, one cannot draw any conclusions from this for the behavior 
of the cut-off errors near Tc. Since the sign of the leading order corrections is negative, the relative 
error quickly develops a pole with growing coupling, beyond which results cannot be extrapolated. 
For longer series the pole might disappear, but the relative cut-off error as a rational function will be 
non-monotonic in general. 



4.9 Comparison of taste splitting for staggered actions 

A main source of systematic error for simulations with staggered fermions is the taste splitting, whose 
effect on thermodynamics has long been underestimated. Recall that the splittings are 

= rriQ + Smf, 5m1 ~ {dsQ^) , (132) 
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Figure 7: Continuum approach of the pressure of massive fermions with a fixed renormahzed quark 
mass mji/T = 0.03. Left: Ideal gas. Right: Two-loop contribution. From [36J . 



to leading order in discretisation errors, where mc corresponds to the true Goldstone boson with taste 
matrix = 75, cf. Table [H In [66] it was observed that the lattice spacing needs to be a < 0.15 fm in 
order for this leading behaviour to dominate, larger lattice spacings still show splittings of even higher 
order in the lattice spacing. A detailed investigation of cut-off dependence of the splitting on finer 
lattices is shown for the asqtad and stout actions in Fig. |H1 with the eight pion multiplets labelled as in 
Tabled] Both actions enter the a^-scaling region within the plot, with the stout action featuring smaller 
splittings and significantly lower mass values for the heavier multiplets on coarser lattices. A similar 
study comparing the scaling between HISQ/tree and stout actions is shown in Fig. [9l with even smaller 
mass splittings for the HISQ/tree action. The lines correspond to linear continuum extrapolations in 
a^a^, where the gauge coupling has been fixed in calculations of the static potential, as = oty [69] . 

Knowing that in the continuum limit the tastes become degenerate, one may identify the particle 
mass with the root mean square average over the tastes, which is another measure for the taste splitting. 
For the example of the pion, 

Mi^*^^ = \\—iMl + M2 + 3M2 + 3M2 + 3M2 + 3M2 + M2 + M^] • (133) 

This quantity is shown in Fig. |9] (right). According to the respective size of taste splittings, this quantity 
is largest for asqtad, whereas the HISQ/tree action features the smallest average mass. Note that the 
differences show up in a region where all actions have an averaged pion mass significantly away from 
the physical value, i.e. the lattices are too coarse to reproduce pion dynamics. In order for the averaged 
mass to be close to the physical mass, lattices with A^,- > 12 are required. 

As discussed in Sec. 13.121 for low temperatures the hadron resonance gas should give a good ap- 
proximation to the equation of state and it is useful to compare its predictions against lattice data. 
However, the taste splitting plots illustrate that for staggered fermions the hadron spectrum is distorted 
by taste splitting and gets fully restored only in the continuum limit. For comparisons with staggered 
fermions at finite lattice spacings, it was therefore suggested [TD] to use the lattice spectrum including 
cut-off effects instead of the continuum spectrum in the hadron resonance gas. For this purpose the 
taste splittings are parametrised in terms of the scale parameter r\ as 

(1 + ci^xr 

An excellent description of the scaling behaviour of the different pion tastes is achieved in Fig. [TOl 
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Figure 8: Pion multiplet masses as functions of the lattice spacing. Left: asqtad action [67]. Right: 
stout action [SH]- Numbers label the taste multiplets, Tabled! The broad band indicates the range of 
lattice spacings for a thermodynamics study at iV,- = 8 between T = 120 and 180 MeV. The narrow 
band in the right panel corresponds to the same temperature range and A^^ = 16. The line labelled "0" 
is the pseudo-Goldstone boson with mass 220 MeV for the asqtad and 135 MeV for the stout results. 
From [681. 




Figure 9: Left: Pion multiplet masses calculated with the HISQ/tree [63] and stout [66] actions as a 
function of a^a^. Right: RMS pion mass with Mq = 140 MeV as a function of the lattice spacing 
for the asqtad ^9], stout [68] and HISQ/tree [63] actions. Vertical arrows indicate the lattice spacing 
corresponding to T ^ 160 MeV for A^,- = 6, 8 and 12. From [63]. 
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Figure 10: The quadratic splittings of non-Goldstone pseudo-scalar mesons in the seven different multi- 
plets calculated with asqtad action at different lattice spacings The lines show the parametrisation 
Eq. (I134p . Open symbols refer to lattice data obtained with the stout action ^66j. From [70] . 

One may now take the non-degenerate tastes of pions and kaons into the hadron resonance gas 
separately with a contribution [70] 

P'"''/T' = ^T7^ E^psln2*^(mpO, (135) 

i=0 

where nips., i = 1 — 7 is calculated according to Eq. f ll32p and riips^ is equal to the pion or kaon mass 
used in the actual lattice calculations. We shall see in Sec. 17.31 that a good description of the asqtad 
equation of state is achieved in this way. 

In summary, taste symmetry violation appears to be a systematic error with a large effect in staggered 
simulations. While formally of 0{asa^) for p4, stout and O^a^a^) for asqtad and HISQ/tree actions, 
its main effect on thermodynamical simulations is in the distortion of the pion spectrum. Identifying 
the Goldstone pion with the physical pion is not sufficient, since also the non-Goldstones contribute to 
the equation of state at finite lattice spacing. Hence, apart from an 0(a^)-errors in thermodynamic 
functions, the taste splitting effectively implies a pion mass out of tune with the physical value and 
distorted contributions of the tastes to the equation of state. 

5 Computing the equation of state by Monte Carlo simulation 

Once a suitable lattice action is chosen, the task at hand is to compute the free energy density from 
the full QCD partition function. A technical obstacle here is that, in a Monte Carlo simulation, one 
cannot compute the partition function directly, since all expectation values are normalised to Z, 

(O) = Z-^Tr(pO) . (136) 

So we need to measure the expectation value of some observable O in order to access the partition 
function. This is achieved by calculating derivatives of the partition function with respect to some 
parameter, as in the derivative and integral methods in the next two sections. 

5.1 The derivative method for anisotropic lattices 

Theoretically the most direct way to obtain thermodynamic functions from lattice simulations is by the 
derivative method using anisotropic lattices with two lattice spacings a^, at in the spatial and temporal 
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directions, respectively. Temperature and spatial volume are then given by 



T = — r^, V = {asNs)\ (137) 



and the anisotropic pure gauge action reads 



S9[U]=Y1 E (l-^ReW^.^x)") 5^ ('l-lReTrf/^4(x)') , (138) 



^s = —,. /3r = ^, (139) 



l<«<i<3 ^ ^ X l<fl<4: 

which is split into terms containing only spatial and temporal plaquettes, respectively. There are now 
two bare lattice gauge couplings, 

2N 2N^ 

oT; Pt ^ 
9s^ 9r 

which are related by the anisotropy parameter ^ = ag/ar- Using the thermodynamic relations Eq. (|2l 
one derives the lattice versions of energy density and pressure as 

Denoting the expectation values of the spatial and temporal plaquettes by 

Rsp =(l- ^ReTiUijix)^ , Rtp=(^l~ ^ReTiUi^ix) ) , (141) 

one gets for the derivatives 

9S\ dPs\^r. f o 9(3. 



{^^ + E + {-^r + a.|^) E Rtp , 

^ ' sp ^ ' tp 



sp ^ ' tp 

Thus, in order to compute the equation of state we need to know the beta-functions of the gauge 
couplings and measure the expectation values of spatial and temporal plaquettes. It is convenient to 
introduce a gauge coupling corresponding to an isotropic lattice with lattice spacing a^. Then the 
anisotropic couplings can be expressed in terms of g and the anisotropy ^, 

gf = g-^ + C,{g',0=9-' + Cs{0 + O{g') 

g^' = g-^ + a{g',0 = 9-' + ct{0 + O{g'), (143) 

where the right hand expressions represent a perturbative expansion for small coupling in lattice per- 
turbation theory [TH [72] . 

Similar expressions can be derived when quarks are included, whereupon also the beta-functions for 
the running of the quark masses appear. Note that also the anisotropy parameter ^ gets renormalised 
and has a corresponding beta-function. These various beta-funcions are only known to some low order 
in perturbation theory, while their accurate non-perturbative determination is complicated. Moreover, 
while anisotropic lattices allow for a finer temperature resolution by choosing small at, there are clearly 
more parameters to tune than in the isotropic case, making for example the determination of lines of 
constant physics more difficult. For this reason, the most popular method in practice is the integral 
method described in the next section and we shall not futher pursue the derivative method here. Recent 
numerical results and a discussion of the associated difficulties using Wilson fermions on anisotropic 
lattices can be found in 
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5.2 The integral method 

The most frequently used method in practice is the integral method [H], in which a derivative of the 
free energy with respect to some parameter serves as observable, which then gets integrated again to 
yield the free energy density, e.g. x = T, 



2^4 



1 

V 



dx 



d_ 

dx 



X 



•\nZ{V,x) 



(144) 



On the lattice, when doing simulations at fixed N^, it is convenient to take derivatives with respect to 
the bare lattice parameters instead, 

{/3,"i/) jyS [■I3,mf 



2^4 



{/3o,m/o) 



/\r3 , 




T=0. 



(145) 



f 

We have subtracted the corresponding vacuum expectation value in each line to renormalise the free 
energy density and pressure to zero for T = 0. The expectation values of the derivatives correspond to 
the average plaquette action and chiral condensate, respectively. 



dlnZ 



N^N^ dl3 
1 d\nZ 




^V/(a;)V'/(5 



(146) 



The vacuum expectation values correspond to A^,- — )■ oo. In a simulation a low temperature lattice 
with sufficiently large A^^ has to be chosen, which introduces an error ~ {Nr/N'^y in the subtracted 
quantity. The integration in Eq. (11451) has to be performed along a line of constant physics, i.e. a curve 
P{mf) such that the renormalised quark masses stay constant. In this case the left hand side corresponds 
to the difference of the renormalised free energy density evaluated at two different temperatures. 



2^4 



{l3o,info) 



f{T) f{To 



-^0 



(147) 



The integration has thus introduced a lower integration constant, which needs to be fixed for the result 
to be meaningful. Note, that for sufficiently small /3o in principle this is possible with an analytic strong 
coupling calculation as in Sec. 13.111 However, these exist only for Wilson's pure gauge action so far, 
whereas dynamical simulations require light fermions as well as modifications for improved actions. 
Many practitioners therefore simply choose /3o low enough so that the pressure is numerically consistent 
with zero and can be neglected. However, this procedure is not practical for the chiral limit, when pions 
become massless and stay relativistic down to very low temperatures. Alternatively, below the quark 
hadron transition the free energy /(Tq) can be matched to the hadron resonance gas, which is again 
validated through the strong coupling expansion, cf. Sec. 13.121 

Finally, using Eq. ([H]), the derivatives appearing in the integrand of Eq. f ll45p are also related to the 
trace anomaly. 



/(T) dT 

2^4 ^ 

m 

rp4 



(148) 
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where the superscript "sub" indicates that the difference with the corresponding vacuum quantity is 
to be taken. We have also converted the temperature derivative into one with respect to the lattice 
spacing, such that we need the lattice beta-functions for the variation of the gauge coupling and the 
masses along a line of constant physics. Hence, the quantity computed directly in a standard approach 
is the trace anomaly, while the pressure and other quantities are obtained by numerical integration of 
those data. Note that this introduces another systematic error, since the trace anomaly data have to 
be interpolated numerically for this purpose. 

5.3 The fixed scale method 

In the last section the temporal lattice extent A^^ was held fixed and temperature was tuned via the 
lattice gauge coupling and the lattice spacing a(/3) in the relation T = {aNr)~^. As a consequence, when 
calculating a given quantity like the thermodynamic functions as a function of temperature, one has 
different lattice spacings and hence different cut-off errors for different temperatures. This commonly 
leads to rather coarse lattices in the low temperature region and to a different quality of continuum 
extrapolation for different temperatures. An alternative use of the integral method is the so-called fixed 
scale approach, where one keeps the lattice gauge coupling, and hence the lattice spacing, constant, and 
varies A^,- instead to tune the temperature. Obviously, this allows only discrete temperature steps and 
the method has only become viable recently, since lattices are now fine enough such that reasonably 
interesting temperature steps can be obtained. 

However, in this case integration over /3 is inapplicable. Instead, one uses T as integration vari- 
able [75]. Starting point is again the evaluation of the trace anomaly followed by integration over 
temperature. 



where b = {P,mf, . . .) is a vector denoting all parameters of the lattice action with the corresponding 
beta-functions db/da along a line of constant physics. The crucial difference is that now integration 
is over temperature while holding /3 fixed, i.e. over A^,-. Clearly, this requires an interpolation to real 
values for A^^, introducing a systematic error which has to be monitored. In [75] it was found that for 
sufficiently fine lattices (corresponding to A^,- > 8 in the usual fixed A^,- approach) results are consistent 
between the two, cf. Sec. O 

The fixed scale approach is complementary to fixed A^^ in that it has large discretisation effects 
at high T, when T > and typical thermal fluctuations occur on the scale of the cut-off. On the 
other hand, at low and intermediate temperatures there are several advantages: for a flxed lattice 
spacing all temperatures are immediately on the same line of constant physics, without any tuning of 
parameters. Furthermore, only one zero temperature simulation at that same /3 is required in order 
to renormalise all observables. Finally, the continuum limit proceeds straightforwardly in complete 
analogy to vacuum calculations, by a series of simulations at different flxed lattice spacings a. Because 
of their complementarity, it will be extremely useful in the future to control systematic errors by using 
them both. 

5.4 Thermodynamic potentials from momentum distributions 

Both, the flxed A^^ and the flxed scale methods discussed above, require the subtraction of vacuum 
expectation values and have to deal with a lower integration constant that is not precisely flxed. These 
circumstances render both methods cumbersome and expensive when either very high or low temper- 
atures are considered. A new approach avoiding these technical difficulties has been proposed in [76] . 




lub 



(149) 
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Starting point is the observation that momentum is a conserved quantity because of translation in- 
variance. Hence quantum mechanical states with different momenta are orthogonal and a projection 
operator P^^^ onto states with momenta p can be defined. The momentum distribution in the grand 
canonical ensemble is then given by 

One may now define a generating function K for the connected cumulants of the momentum distribution, 



-l^e^P-i?(/3,/x,p), 
p 

K{2nu2nM 0^^^ d^^^ d^^^ K{z) 



P 



(_l)ni+n2+n3 + l ^Z^^^ (9z^"^ ^zf ^ 



(151) 

z=0 



Employing Ward identities for the conservation of energy and momentum, one can show that for fi = 
the entropy s = {e + p)/T is connected to the second order cumulant, while the specific heat is related 
to a combination of the second and fourth order cumulants [761 EZj, 

^2,o,o((T, /i) = T(e + p) , = = -jr^ . (152) 

On the other hand, the generating functional for the momentum distribtuion can also be expressed as 
a ratio of partition functions, 

Z (T, /i) 

where Z{T,ii,z) is a partition function in which states of momentum p are weighted by a phase e*P'^. 
It can be represented as a Euclidean path integral, where the weight factor is absorbed in the shifted 
boundary conditions for the generic field, 

0(1/T,x) = ±0(O,x + z) . (154) 

Note that the cumulants represent connected correlation functions of conserved charges, i.e. momentum. 
As such they are physical observables and do not require renormalisation, even on the lattice [77], 
which therefore saves the vacuum subtraction required in the integral method. Hence, calculation of 
the ratio of the shifted and unshifted partition functions permits an evaluation of the entropy and other 
thermodynamic functions of relativistic field theories. Applications to SU{3) pure gauge theory and 
scalar field theory have been presented in [761 EZ] , respectively. 

For pure gauge theory, the required ratio of partition functions is split up in factors 

Z{T) -llz(T,r.^O' ^ ^ 

with partition functions Z{T, rj), (r^ = z/n, i = 0, 1, . . . , n) and corresponding actions 

5(f/, n) = nS{U) + (1 - r,)S{U') , (156) 

where obeys shifted boundary conditions, see Eq. fll54p . The continuum value of the second cumulant 
is obtained from a discrete lattice by a limiting procedure, 

s{T) _ K,MT) _ 2K{T,z,a) 
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where z = {uza, 0, 0) and Uz has to stay constant when a — )■ 0. Numerical results obtained by this 
method are shown in Sec. 16.21 Note that, while no subtraction of vacuum lattices is necessary in this 
approach, it requires the computation of an entire series of products of partition functions, Eq. (I155p . 
which is again computationally costly. It will be interesting to see in future investigations whether 
an optimised choice for the interpolating partition functions can render this method computationally 
advantageous. 



6 Numerical results for pure gauge theory 
6.1 SU{3) Yang-Mills 

The results for the equation of state of pure gauge theory [78] are by now "historical" and well-known. 
Since systematic errors in this case are relatively benign, it still represents the only case with a complete 
continuum extrapolation and sets the standard of what can be achieved with lattice simulations. It 
is based on Wilson's standard pure gauge action using lattices with A^^ = 4, 6, 8 and aspect ratios 
NjNr = 4 - 5.3. 

SU{3) pure gauge theory has a first order finite temperature deconfinement transition at unique 
critical couplings f3c{Nr) associated with centre symmetry breaking. Sec. 12.21 In [78] those were deter- 
mined from the location of the peak of the Polyakov loop susceptibility and converted to units of the 
string tension by calculating the latter on 32^ vacuum lattices for the same /5- values. In the continuum 
limit, 

T 

= 0.625 ± 0.003 (+0.004) , (158) 

^/a 

where the second error corresponds to estimated finite size corrections to the critical couplings. Together 
with s/a = 420 MeV this gives a critical temperature f» 265 MeV in terms of which all other 
dimensionful quantities can be expressed. 

The equation of state was calculated by the integral method. Sec. 15.21 The trace anomaly for the 
pure gauge case reduces to 

^ = N'a^ (3(Trf/; + Trf/^) - Q{TrU,)T=o) , (159) 

and requires evaluation of temporal and spatial plaquettes as well as the lattice beta-function. For the 
latter the difference in critical couplings corresponding to a doubling of the lattice spacing is computed, 
= /3c(2A^r) — (3c{Nt-), and modeled by the perturbative two-loop beta-function modified by an 
effective coupling and a scaling factor. The corresponding trace anomaly for the three different values 
of Nj. is shown in Fig. [11] (left), including smooth spline functions as numerical interpolation between 
the data points. Note that, for a first order phase transition, a discontinuity in the energy density is 
expected in the thermodynamic limit. However, on a finite size lattice, there is tunneling between the 
metastable phases resulting in a continuous average value for all temperatures. The pressure is then 
obtained by integrating the trace anomaly according to the first line of Eq. fll45p . The free energy 
density is neglected (statistically compatible with zero) for a lower integration limit below the critical 
coupling, /3o < Pc- 

The numerical results show significant cut-off effects for = 4. A continuum extrapolation is thus 
performed for the pressure using the A^^ = 6, 8 data based on the ansatz 

with a fit parameter c{T). The resulting continuum curve, as well as those for energy and entropy 
densities, are shown in Fig. [11] (right). 
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Figure 11: Left: The trace anomaly (e — 3p)/T^ for pure gauge theory calculated on lattices with 
temporal extent Nr = 4, 6, 8. Right: Continuum extrapolated results for pressure, energy and entropy 
density. From [78] . 



These results have proven stable against checks with improved actions, Fig. [12] (left). The data show 
that both the tree-level improved and tadpole improved Wilson actions used in [79j approach continuum 
results for the pressure already on rather coarse lattices with A^,- = 3, 4. In Fig. [12] (right) the trace 
anomaly is compared between the A^,- = 8 data from [78] and those produced with the fixed scale 
approach in [75]. In this work, the standard Wilson action was used with fixed lattice gauge couplings 
P = 6.0,6.1. Using tq to set the scale, this corresponds to the lattice spacings a = 0.093,0.068, 
respectively. Temperature was then tuned by varying N^- leading to the data points in the figure. 
For the zero temperature subtraction, A^,- = 16,22 lattices were used for the i2,i3 runs in the figure, 
respectively. For the beta-function, the data from [78] were used. 

Fig. [12] (right) shows generally good agreement between the two methods. At low temperatures, 
one observes a finite volume effect, but no difference between the lattice spacings. On the other hand, 
at large temperatures, the fixed-scale results overshoot the fixed A^^ results. This is likely because for 
the former the cut-off effects are stronger at high temperatures, as discussed in Sec. 15.31 A sensitive 
measure for cut-off effects is provided by the peak height of the trace anomaly. It indicates that the 
continuum limit has indeed been reached in Yang-Mills theory, whereas for full QCD the peak height 
is not yet settled between different actions, as we shall see later. 



6.2 5*^7(3) at high temperatures 

As discussed in Sec. 13. 8[ simulations of thermal systems require 1 <^ A^^ <^ A^s in order to have man- 
ageable cut-off and finite size effects. This has limited simulations of the QCD equation of state to 
typically T < 5Tc. On the other hand, one would like to make contact to weak coupling perturbation 
theory, which unfortunately converges rather poorly until much higher temperatures are reached. Im- 
provements in the convergence behaviour can be achieved by employing screened perturbation theory, 
consisting of partial resummations to infinite orders rather than strict perturbation theory, cf. [80]. Re- 
cent calculations using the HTL-scheme have managed reasonable contact with lattice results at ~ 4Tc 
for pure gauge theory [81] and dynamical QCD [82]. Nevertheless, it still is desirable to extend lattice 
simulations to higher temperatures where perturbation theory is fully controlled. 

One possibility to achieve this is by use of the momentum distribution method discussed in Sec. 15.41 
Its feasibility for pure gauge theory has been demonstrated in [76], with numerical results obtained 
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Figure 12: Left: Comparison of pure gauge calculations with tree-level, tadpole improved and a fixed 
point action on Nr = 3 and 4 lattices. The arrows indicate the ideal gas result in the continuum limit. 
The solid line is the continuum extrapolation using the unimproved Wilson action [78]. From [79] . 
Right: Trace anomaly from the fixed scale approach [75], compared to the = 8 data from [78] . The 
dotted lines are spline interpolations, horizontal errors from insecurities in the scale are smaller than 
the symbols. From [75] . 



from the standard Wilson action as shown in Fig. [12] (left). Simulations have been performed for three 
different temperatures and the scale was set using ro and its related beta-function from |83j . Within the 
uncertainties good agreement for the lower temperatures with earlier results from the integral method 
[78] is observed. The continuum value for the entropy has again leading O(a^) corrections. 

Another attempt to develop methods in this direction has been reported in [HU |85], based on a 
modification of the integral method. Starting point is the observation that the vacuum divergence, 
which is subtracted to renormalise the free energy density or pressure, is temperature independent. 
The subtraction is therefore also achieved at two non- vanishing temperatures. In [HI] it is suggested to 
consider the series with half-temperature steps, 

p{T)ren ^ p{T) - p{0) _ p{T) - p{T/2) + p{T/2) - p{T/A) + p{T/A) ... - p{0) 

p{T)-p{T/2) 1 p(T/2)-p(T/4) 1 p(T/4) - p(T/8) 

^ 16 (T/2)4 ^ 256 (T/4)4 ^ ' ' ' ^ ^ 

The idea is to evaluate the different terms using only two lattices with and 2N^. The temperature 
for the successive terms is lowered by tuning the lattice gauge coupling. In practice, only a few terms 
are needed, because of the rapid suppression from the prefactor. This is demonstrated in Fig. [T3] (right), 
which is based on the data from = 6, 8 lattices [78], i.e. with a suppression factor 6/8 instead of 2. 

However, a word of caution is in order. Eq. (I16ip is written in continuum notation, where each 
term is represented by a /3-integral from the lower integration limit /3o to some /3(T/x, A^,-). The upper 
integration limits, and hence lattice spacings, are different in each term, leading to different cut-off 
effects. Thus, at finite lattice spacing the added and subtracted terms do not cancel exactly, as the 
continuum formula suggests. Rather the continuum limit should be taken for each term separately. It 
would be interesting to further investigate how improvement alleviates this problem or how small lattice 
spacings need to be to neglect this effect. 

This proposal was applied to a calculation of the equation of state for pure gauge theory at high 
temperatures using a tree- level Symanzik improved action [85]. The savings on large lattices were 
then invested in a large range of aspect ratios Ng/Nr = 4 — 24 and dedicated finite size analyses. Fig. 
(left) shows the trace anomaly calculated on lattices with Nr = 5, 6, 7, 8, in comparison with the known 
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Figure 13: Left: Scaling behavior of s/T^, Eq. f ll57p . The Stefan-Boltzmann value reached in the high-T 
limit is also displayed. From [TH]. Right: Pressure in lattice units calculated on 32^ x 6 and 32^ [75] 
compared to the sum of individual terms "p' of Eq. f ll6ip using 32'^ x 6 and 32^ x 8 lattices. From |84] . 



low temperature lattice data from [78] as well as with standard perturbation theory through order [86] 
and HTL-resummed perturbation theory [81]. Note that the dimensionless trace anomaly I{T)/T^ has 
been rescaled by (T/T^)^, in order to show deviations from perturbative behaviour. Since the effective 
expansion parameter for the latter is g'^T/niE, where ~ gT is the electric screening mass scale, at 
most logarithmic corrections in temperature can modify the behaviour of the dimensionless trace 
anomaly in perturbation theory. Hence the linear, nearly constant section up to ~ 4Tc is indicative of 
non-perturbative effects in the pressure, which then weakens at higher temperatures. 

A study of the volume dependence and hence the requirements on the aspect ratio was also performed 
on A^T- = 5 lattices. From Ng/Nr = 6 upwards the data do not show any finite size effects. This is an 
important observation, since the soft modes of the magnetic gauge fields live on the scale ~ g'^{T)T and 
for high temperatures correspond to a large correlation length l/{g'^{T)T). However, the corresponding 
screening masses are those of the lightest glueball. However, the coefficient is large, mo++ ~ 2.5g^{T)T 
[57] , ^ ~ OA/g'^{T)T and T grows faster than g{T) shrinks. It is thus plausible that the finite size 
requirements are fulfilled on modest lattice volumes as observed in the simulations. 

A continuum limit was then taken for the aspect ratio Ng/N^- = 8 and the corresponding pressure is 
shown in Fig. [H] (right). The perturbative regime can thus be reached at about T > lOTc and the data 
can be used to fit the (/^-coefficient in the perturbative expression for the pressure p6]. The dashed 
curves in the figure respresent a modelling of the non-perturbative temperature regime Tc < T < 4Tc, 
by adding to the perturbative formula up to g^ the contributions, 

P = Ppert(/)+anp , ttnp = 0.879(2) (40) 

P = Ppert(/) + anp + fcnpln(T/T,) a,p = 1.371(1)(50), ft^p = -0.618(2)(4) . (162) 

The parametrisation including the log gives a good description of the data all the way down to the critical 
temperature. Together with the perturbative expression from [86], we then have a parametrisation of 
the equation of state for all of the deconfined phase at our disposal. Note that even at temperatures 
T ~ lOTc the pressure is still ~ 10% below the ideal gas limit. This is of course consistent with the 
logarithmic running of the coupling constant. It illustrates however, that there are still significant 
interactions of the soft modes in the plasma and one must be very careful with intuition based on the 
Stefan-Boltzmann limit. The latter becomes accurate only at exceedingly high temperatures. 
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Figure 14: Left: Trace anomaly for different compared with previous lattice results [78], perturbation 
theory [86j and HTL-resummation [81]. Right: Continuum extrapolated pressure normalised to the ideal 
gas value and comparison with perturbative expressions. From [85] . 




SU{N) pure gauge theory is also of theoretical interest in the limit — )■ oo. Firstly, one would like to 
understand how far SU{3) is from this limit, for which considerable calculational simplifications apply 
and analytic treatments become feasible [88j- Secondly, the large N limit is a crucial ingredient for the 
modelling of QCD-like theories based on the AdS/CFT conjecture p9]. In order to see to which extent 
SU (3) can be described by large theories, a non-perturbative calculation of the equation of state by 
lattice methods is useful also for > 3. 

Calculations for the critical temperature Tc have been extended to A^ = 10 [9Qj, the equation of 
state exists for A^ = 3,4,5,6,8 [9ll [92]. The simulations for the latter were run on 20^ x 5 lattices 
for A^ = 3 and 16^ x 5 lattices for A^ > 3. The calculation employs the integral method and the scale 
for SU (3) was set using tq [83] . For larger A^, it was based on interpolated string tension values from 
[93] . The numerical results are displayed in Fig. [151 A rather remarkable scaling is observed for all 
thermodynamic potentials. Also shown are results from a prediction from a holographic QCD model 
[M] which agrees with the data remarkably well. 

7 Numerical results with staggered fermions 

As described in Sec. l3.5[ staggered fermions are computationally the cheapest because they correspond to 
one-component fields in the action Eq. ( I44p rather than four-spinors. Their leading discretisation errors 
in the equation of state are C(a^) in the standard formulation and suppressed when using improved 
actions. For those reasons simulations using staggered fermions are most advanced, in the sense of 
being closest to the physical masses and using the finest lattices. Extended series of simulations have 
provided evidence that, at the long last, systematic errors are getting under control and which actions 
have the best scaling behaviour in a given temperature range. 

7.1 The pseudo-critical temperature 

Let us recall that, for QCD with physical masses and a large neighbouring region, there is no chiral 
or deconfinement phase transition, merely an analytic crossover. Hence there is no uniquely defined 
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Figure 15: Trace anomaly and pressure, obtained for various numbers of colours, N . From |92l 



critical temperature, only a pseudo-critical "Tc" depending on the observable used to define it. 

The determination of the pseudo-critical temperature for QCD at the physical point has been the 
subject of some apparently contradictory results. Based on simulations with the asqtad action on 
A^^ = 4, 6, 8, T, = 169(12) (4) MeV was quoted [95]. This was followed by the prediction = 192(7) (4) 
MeV based on the p4 action on A^^ = 4,6 [HS] • In both cases simulations were performed at several 
light quark masses larger than physical and then extrapolated simultaneously to the physical point and 
the continuum limit. In [96] the extrapolation formula is 

r,T,{mu mf y^ N,) = rMT^t''^ oo) + A{r^m^,Y + ^ . (163) 

T 

Here 0(4) scaling in the chiral limit is assumed, the data being consistent with the associated crit- 
ical exponent d. Moreover, the values for Tc based on susceptibilities of the Polyakov loop or chiral 
observables were found to be within 185 — 195 MeV, the quoted value representing an average. 

By contrast, in a series of simulations using a stout link improved staggered action over a wide 
range of = 4, 6, 8, 10, 12, 16 [211 [66l [68|, consistently lower values were quoted. The simulations were 
performed on the physical point and hence no extrapolation in the quark mass was necessary. For the 
continuum extrapolation the four largest A^r corresponding to the finest lattices were used, providing 
moreover different values when extracted from the chiral condensate A; Eq. (162|) . or the Polyakov 
loop, Eq. dSD, [SHlEn]: 

r 157(3) (3)MeV A,,. 
' \ 170(4) (3)MeV L ' ^ ' 

The continuum extrapolation is based on four lattice spacings, which allows to monitor that cut-off 
effects are small and truly scale as ~ 1/A^r- Note that the observable-dependent difference in is 
completely consistent with the fact that there is no true phase transition at all, neither a chiral nor 
a deconfining one, merely a transition region in which the chiral and confinement properties change 
gradually and smoothly. 

In the meantime the differences between these actions are understood and nearly resolved, especially 
since a new set of simulations with asqtad and HISQ/tree actions on A^^ = 8, 12 and an additional lighter 
quark mass [63] is available. A combined chiral and continuum extrapolation, based on the observable 
Ai^s gives Tc = 154(9) MeV, in complete agreement with the stout results above. The previous differences 
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are largely due to underestimated cut-off effects due to taste splitting in the asqtad and p4 actions. 
Moreover, the corresponding data obtained for = 4, 6 are not in the A^~^ scaling region and could 
not yet be extrapolated with the formula used. The same effects are visible in the equation of state, and 
we shall discuss them together when we compare results between different actions in detail in Sec. [91 

7.2 Flavour dependence and interpretation of the equation of state 

The results of a computation of the pressure with the p4 action [HZ] are shown in Fig. (left). The 
data have been obtained for Nf = 2,3 with (bare) mass rUq/T = 0.4 as well as for Nf = 2 + 1 with 
a heavier mass rriq/T = 1 on a very coarse lattice, Nj- = 4. For comparison, continuum extrapolated 
pure gauge results are also included. Without entering the details of the calculation, the plot tells us 
some interesting qualitative features. As in the pure gauge case we see a rapid rise of the pressure in a 
narrow transition region. The pseudo-critical temperature as well as the magnitude of p/T^ rise with 
the number of degrees of freedom liberated at the transition. This last conclusion is firm, since the 
pressure also rises for fixed temperature when light quarks are added to the theory, consistent with the 
expectations based on the Stefan-Boltzmann limit. We may thus conclude that the rapid rise in the 
pressure signals deconfinement, in the sense that there are more and lighter degrees of freedom in the 
system. On the other hand, these are not yet free quarks and gluons, as the pressure falls significantly 
short of its Stefan-Boltzmann limit. 



7.3 Equation of state for almost physical quark masses: hotQCD 

The hotQCD collaboration employs and continues previous work of the MILC and RBC-Riken collab- 
orations. They have been working with asqtad and p4 actions over several years. Simulations exist for 
Nf = 2 + 1 using three different lattice spacings, with A^^ = 4, 6 published in [98] (asqtad) and [28] 
(p4), while combined results for both actions and = S have been presented in [TT] . 

The actions used for those simulations are those specified in Sec. 14.61 The strange quark mass rris is 
fixed close to its physical value while ms/mud = 10. This corresponds to a pion mass of m^r ~ 220 MeV. 
The lines of constant physics are defined by keeping the meson masses in Table H] fixed, where the scale 
is set by the Sommer parameter tq. The latter is related to the T(2S' — IS*) splitting in calculations 
with the asqtad action [^3 1100] . This gives tq = 0.469(7) fm. For = 8, an aspect ratio Ng/N^- = 4 
was used throughout. 

The beta function for the lattice gauge coupling was fixed by fitting the /3-dependence of vq to the 





p4-action 


asqtad action 


mss/rriK 


1.59(5) 
1.33(1) 
0.435(2) 


1.83(6) 
1.33(2) 
0.437(3) 


m^r-o 


0.371(3) 
1.578(7) 
1.158(5) 





Table 4: The strange pseudo-scalar mass = ^J2m?j^ — in units of tq and mass ratios that 
characterize lines of constant physics for fixed ratios of light and strange quark masses, mi/ms = 0.1, 
above, and mi/ms = 0.05, below. Errors represent the variation over the lines of constant physics. 
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Figure 16: Left: Flavour dependence of the pressure for A^.^ = 4 lattices compared to the continuum 
pure gauge result. From [97] . Right: The trace anomaly for mi/ms = 0.1 calculated on lattices with 
temporal extent A^.^ = 6, 8. From |17] . 



ansatz 

ro ^ 1 + era\/3) + fra\/3) 

a a,i?2(/3) (1 + bra\(3) + c,a4(/3) + rf,a6(/3)) ' ^ ^ 

where R2{P) is the two-loop beta function for 3- flavor QCD and a{(3) = i?2(/3)/-R2(3.4). Similarly the 
beta function for the quark mass was obtained from a parametrisation of the bare quark mass 

miro/a = m^'^'ro i^) P(/3) , (166) 

with a sixth order rational fit function P{/3) and the renormalisation group invariant m^'^^ . 

The trace anomaly calculated for both types of actions is shown in Fig. Hn] (right). Good agreement 
between different lattice spacings and actions is observed for low and high temperatures, but differences 
remain in the peak region, indicating that there are still significant discretisation effects. In order to 
obtain a smooth curve, the data are interpolated numerically. The remaining thermodynamic quantities 
are then obtained by numerical integration, Eq. f ll45p . and using the thermodynamic relations Eq. (|2]). 
This gives the results for the pressure, energy and the entropy. Fig. [TTl and the speed of sound. Fig. US] 
(left). The integration requires fixing the lower integration constant by normalising the pressure at 
some low temperature to a prescribed value. The figures correspond to setting p = at Tq = 
and integrating the fits to the trace anomaly. Normalising the pressure to the hadron resonance gas 
prediction at Tq = 100 MeV instead leads to a constant shift marked by the black bar in the upper 
right of the figures. The shape of the velocity of sound curve appears to depend on the details of the 
numerical interpolation of the data for the trace anomaly, as is evinced by the difference between the 
data points and the smooth curves. This difference should be interpreted as a measure of the systematic 
error associated with the processing of the original data for the trace anomaly. 

These results were supplemented subsequently with an additional set of simulations using the p4 
action with lighter quark masses, mi/nis = 0.05, corresponding to m^r ~ 156 MeV, using 32^ x 8 finite 
temperature and 32"^ zero temperature lattices |101] . The line of constant physics in this case was 
defined by the second set of quantities in Table HI The resulting trace anomaly is shown in Fig. [TS] 
(right) and compared to the previous calculations for heavier light quark masses. In order to see the 
mass effect, the Nr = 8 data should be compared. Differences are visible for T<200 MeV, where the 
lighter quark values sit systematically higher. This is consistent with the decrease of the pseudo-critical 
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Figure 17: Left: Asqtad and p4 results for nii/ms = 0.1 on A^,- = 8. Crosses with error bars indicate 
the systematic error from different integration schemes. The black bar at high temperature indicates 
the systematic shift of data when matching to a hadron resonance gas at T = 100 MeV. The vertical 
band indicates the transition region 185 MeV < T < 195 MeV. Right: Entropy density. From |47] . 




Figure 18: Left: Pressure divided by energy density (p/e) and velocity of sound (c^) Lines without 
data points give calculated from the interpolating curves for e/T^ and p/T^. From [17]. Right: 
Comparison of mi = O.lm^, 0.05ms as well as with the HRG model which includes all the resonances 
up to 2.5GeV. Also shown are interpolations of the lattice data. From |101] . 
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Figure 19: Left: Meson masses and ratios of quark masses and decay constants on the finest lattices 
denoted by squares in the left plot. Horizontal bands give experimental results. From [66j. Right: Lat- 
tice spacing as a function of /3 for A^j = 2 + 1 simulation. The curve corresponds to 2-loop perturbative 
running. From |1U3] . 



temperature for the crossover (within a given definition) with decreasing quark mass. On the other 
hand, for temperatures above the transition the quark mass dependence largely disappears. This is 
consistent with the fact that hadrons are deconfined and the screening masses corresponding to the 
same quantum number channels are proportional to ~ 27rT, compared to which the bare quark masses 
are negligible. 

Note that, for T < 180 MeV, data from both masses are significantly lower than the hadron resonance 
gas prediction, pointing at the presence of discretisation errors in the light hadron sector. This is 
corroborated by a comparison with hadron resonance gas predictions using the distorted lattice spectrum 
including cut-off effects such as taste splitting [70], as outlined in Sec. 14. 9[ and is shown in Fig. |29] (right) 
in Sec. [H When the hadron resonance gas model includes the taste splitting and cut-off effects in the 
mass values, excellent agreement is found between its predictions and the lattice data for T < T^. One is 
thus led to conclude that the observed deviations from the continuum hadron resonance gas predictions 
is to a smaller extent caused by larger than physical quark masses, and to a large extent by discretisation 
effects in the low lying hadron spectrum. In the final limits of physical masses and vanishing lattice 
spacing one would thus expect agreement with the continuum hadron resonance gas predictions. 



7.4 The equation of state for physical quark masses: BMW 

At the time of writing, only the Budapest-Marseille- Wuppertal collaboration has provided a series of 
calculations for physical quark masses. They are also working on the finest lattices used in this context 
and provide a full continuum extrapolation for a growing number of temperature values. Moreover, 
comparison between different lattice spacings suggests that cut-off effects are indeed small. Results for 
the equation of state on A^^ = 4, 6 can be found in |102] . A^^ = 8, 10, 12 were added in |103[ 1105] . 

The actions used in this series of simulations are a tree-level Symazik improved gauge action. Sec. 14. 2[ 
and the one-link staggered action, Eq. (jH]), where the link is stout-smeared. Sec. 14.61 The link smearing 
was applied twice, with parameter p = 0.15. The bare lattice parameters are fixed using as input the 
physical values m.,^ = 135 MeV, fx = 155.5 MeV, rriK = 495 MeV [TT] to define the line of constant 
physics. This translates to a fixed bare quark mass ratio ms/rriud = 28.15(1) along the line ms{P). 
The required beta functions were computed in two different ways. For 3.45 < (3 < 3.85 the hadronic 
observables were calculated on zero temperature lattices 24^ x 32 up to 48^ x 64, setting the scale using 
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fx- The masses of the fl, 0, K* mesons were then monitored to stay at their physical values, Fig. [19] 
(left). For the larger couphngs 3.9 < (3 < 4.6 corresponding to the smallest lattice spacings, spectrum 
calculations would have been too expensive. Instead Creutz ratios were used to define an effective pure 
gauge coupling [104j . which was then modified by including perturbative one- loop running of the quark 
masses. For a few f3 > 3.99, control simulations were run to ascertain that the deviation from the true 
line of constant physics in this way remained smaller than 2%. The lattice gauge couplings employed 
reach the weak coupling region, where the corresponding beta-function is described by perturbative 
two- loop running. Fig. [19] (right). 

A check for finite volume effects was run for A^,- = 6 with aspect ratios Ng/Nj. = 3, 6, corresponding 
to box sizes of ~ 3.5,7 fm, respectively. In addition, aspect ratios Ng/Nj. = 3,4 were compared on 
Nt- = 8 [lU5j. No differences beyond the statistical errors were found, so all simulations on finer lattices 
were performed with Ng/Nr = 3. 

The main results characterising the equation of state are shown in Fig. [201 where tree-level improve- 
ment on the level of the observable has been applied in order to reduce cut-off effects. For the trace 
anomaly, there are data with four different lattice spacings. The = 8, 10, 12 data essentially lie on 
top of each other, indicating that cut-off effects are very small in this regime of lattice spacings. Only 
the Nt- = 6 data show some deviation on the right flank of the peak. The inset shows a comparison 
with the continuum hadron resonance gas model using physical masses, which describes the data well 
up to T ~ 140 MeV. The panel showing the pressure consists of three lattice spacings corresponding to 
Nr = 8, 10, 12. The two coarser lattice spacings cover a temperature range over an order of magnitutde, 
T ~ 100 — 1000 MeV. Since the lattice spacings shrink with T, the results are most reliable at high 
temperatures. Interestingly, A^,- = 6,8 show no difference at higher temperatures and one may consider 
these lattices fine enough, even though stout fermion actions are not optimised for the ideal gas. This 
is consistent with the observation that for the highest temperatures, the pressure is still significantly 
short of the Stefan-Boltzmann limit. This means interactions still play an important role for those tem- 
peratures, even though the rise of the pressure signals deconfinement. Of course, this is to be expected 
based on the slow logarithmic running of the coupling as{T). 

On the other hand, at the lowest temperature T ^ 100 MeV, the lattice spacing is a ~ 0.2 fm on 
At = 10, which is too coarse to be in the scaling region. This is also indicated by the fact that for this 
lowest temperature the pressure from the lattice is only about two thirds of what the hadron resonance 
gas predicts. While the latter of course does not constitute a complete result either, it should be better 
the lower the temperature. The authors quote the discrepancy to the hadron resonance gas as a measure 
of the remaining systematic error for low temperatures. Simulations at an additional At = 12 followed 
by a continuum extrapolation should help to provide a final answer, see below. Finally, Fig. [20] also 
shows the entropy and energy densities, while the speed of sound is depicted in Fig. [21] (left). 
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6.79 
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5.59 


7.34 


-5.60 


1.42 


0.50 



Table 5: Parameters of the function in Equation ( I167p describing the trace anomaly in the Nf = 2 + 1 
and in the Nf = 2 + 1 + 1 flavor cases. 

Note that the results presented in the plots do not yet constitute a continuum extrapolation over the 
entire temperature range, which are however provided for T = 132, 167, 223 MeV. Nevertheless, on the 
grounds that no systematic shift with lattice spacing is visible, the authors provide an estimate for the 
continuum result by averaging the two finest lattices A^t = 8, 10, with the low temperature discrepancy 
from the hadron resonance gas result quoted as a systematic error over the whole range. The resulting 
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Figure 20: Left: Equation of state as a function of temperature on = 6, 8, 10, 12. Upper Left: Trace 
anomaly. Upper Right: Pressure. Lower Left: Energy density. Lower Right: Entropy density. From 



trace anomaly is shown in Fig. [21] (right), together with a fit function which provides a convenient 
closed-form parametrisation over the entire temperature range, 

HT) ( u u u /,2^ (. , /o-[tanh(A-t + /2) + l] ^ 

Here, t = T/(200MeV) and the parameter values describing the simulation data are given in Table [51 

Since the publication of [103] . further data were collected on Nr = 12 lattices, working towards a 
true continuum extrapolation, with preliminary results presented in [105] . The data points for all A^r 
were interpolated by standard cubic splines and then extrapolated to the continuum using a ~ l/iV^ fit 
ansatz. The results are shown in Fig. [22] (left). Note that tree- level lattice corrections were subtracted 
from the plotted data, but not from the data entering the continuum extrapolation. Excellent agreement 
with the hadron resonance gas prediction using physical masses is obtained. 

Finally, for purposes of comparing with simulations at different parameter values also the effect of 
quark mass has been studied on an A^,- = 8 lattice. This is shown in Fig. [25] (right), where the results for 
the physical parameter values for all three quarks are compared with the mass-degenerate Nf = 3 case, 
where ,i = mf^^^. This corresponds to a rather heavy m-Tr ~ 720 MeV. Qualitatively, similar changes 
as in Fig. [18] are observed, only much more pronounced as the mass difference is larger. The peak shifts 
to lower temperatures and decreases in height with shrinking light quark masses. Note that for large 
temperatures the quark mass dependence disappears. This is consistent with the fact that confinement 
is lost while the leading order screening masses in the corresponding quantum number channels are 
~ 27rT, against which the bare quark masses are negligible. This is also why it makes sense to compare 
with massless perturbation theory in the plot. 
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Figure 21: Left: Speed of sound on A^^ = 6,8,10. Right: Trace anomaly based on the average of 
Nr = 6,8 data as "continuum estimate". Also shown is the parametrisation Eq. f ll67p . From [103] . 




Figure 22: Left: Attempt for a continuum extrapolation based on four lattice spacings. The plotted 
data points are tree-level improved, but not the data entering the continuum extrapolation. From [105] . 
Right: Mass dependence of the trace anomaly. From [103] . 



7.5 Towards inclusion of the charm quark 

Several years ago, based on next-to-leading order perturbative investigations it has been suggested that 
the equation of state gets modified by thermalised charm quarks at temperatures starting at T > 2Tc 
|106] . as depicted in Fig. [22] (left). If thermalisation happens sufficiently fast, this would affect heavy ion 
collisions at LHC, where such temperatures will be reached, and certainly the equation of state of the 
Standard Model in the early universe [2] . The behaviour predicted in perturbation theory was confirmed 
by initial lattice studies where the charm was treated in the quenched approximation [107[ 1103] , Fig. |23] 
(right). One would imagine this to be a good approximation, since quantum fluctuations of the charm 
(quark loops) are suppressed by its mass. On the other hand, the same argument would suggest a small 
effect on the physics, so the question merits a closer look. 

A detailed investigation, comparing quenched and dynamical simulations, was recently started in 
[105] . According to the formula used for the integral method, Eq. fll45p . the charm quark may modify 
the expectation values of the gauge action and the chiral condensate. Its influence on the beta-function 
is known to be small. Fig. [2l] (left) compares the contribution of the gauge action with and without 
a dynamical charm quark, with practically no difference visible, and the same flnding for the quark 
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Figure 23: Left: Perturbative calculation of the pressure matched to the resonance gas for three and 
four flavours of quarks. From |106j . Right: Lattice results for Nf = 2 + 1 + 1 with charm in the quenched 
approximation compared to Nf = 2 + 1 (dynamical) |1U7] . From |1U8] . 
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Figure 24: Left: Gauge action used in the integral technique from the Nf = 2 + 1 data set, as in 
a partially quenched calculation, compared to the preliminary results for Nf = 2 + 1 + 1. Right: 
Nf = 2 + 1 + 1. The data points are tree- level corrected and compared to the Nf = 2 + 1 continuum 
estimate. From |1U5] . 



contribution to the pressure. The effect of the charm quark on the pressure then enters likely through 
the integration, i.e. solely through the modification of the line of constant physics. Nevertheless, in 
|1U5] the full dynamical calculation was performed, with its result in Fig. [21] (right). 

The result shows clear differences to the quenched calculation discussed above. The main observation 
is that the effect of the charm quark is delayed to T ~ 300 MeV now, similar to what is seen in the 
perturbative estimate. This is a very interesting observation, calling for a closer investigation. In 
particular, one would now like to see an explicit evaluation of the effect of the charm on the lines of 
constant physics. Another interesting question concerns cut-off effects. It is well known that at zero 
temperature heavy quarks are beset by large discretisation errors when their mass gets too close to the 
cut-off. By contrast, the data in Fig. |2l] (right) appear to show no significant cut-off effects even at low 
temperatures, where the lattices are coarsest. This is somewhat surprising, since a A^^ = 6 lattice at 
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low temperatures corresponds to a>0.2 fm. At high temperature the lattice spacing is finer and the 
situation should be better. It will be interesting to see future investigations of these questions. In any 
case, the conclusion emerges that for temperatures T>400 MeV indeed a thermalised charm quark 
significantly affects the equation of state of QCD and the Standard Model. 



8 Numerical results with Wilson fermions 

The Wilson discretisation of the fermionic action gets rid of the doublers dynamically, by providing them 
with masses ~ which lead to decoupling as the continuum is approached. There is no reduction 
of spin degrees of freedom on the level of the partition function, as for staggered fermions, and hence 
the simulations require both more memory and time. To a large part this explains why there are less 
Wilson results at finite temperature. For a long time, a study of the equation of state existed only for 
two fiavours on relatively coarse lattices |109] . Recently some efforts were made to study the quark 
hadron transition on finer lattices up to A^^ = 12 and lighter quark masses, down to rrij^ ~ 400 MeV 
|110[ IllH 1112] . but these are again for Nf = 2 and so far without results for the equation of state. The 
first thermodynamics studies with Nf = 2 + 1 fiavours of Wilson fermions appeared only very recently 
[TT31I261I27]. 

In jll3j the trace anomaly was calculated based on the equation, 

e-3p Nf(dl3/dS\ , dK^a / dS \ , dK^/dSX \ 



T'^ V da\d(5 / da \ dn^d /sub \ 



sub 



The authors work with the fixed-scale approach, cf. Sec. 15. 3[ using an RG-improved gauge action. 
Sec. 14.21 and non-perturbatively clover-improved Wilson fermions. Sec. 14.31 For the finite temperature 
configurations, lattice 32^ x Nr with A,- = 4, 6, 8, . . . 16 are used. For the zero temperature subtraction, 
they use data from the CP-PACS+JLQCD collaboration [114] . These correspond to a physical strange 
quark mass, defined by m^^^/m$ ~ 0.74 and light quark masses defined by m^/mp = 0.63, with pions of 
~ 620 MeV. On the 28^ x 56 vacuum lattices, the lattice spacing was determined as ro/a = 7.06(3) 
[TT5] . or a ^ 0.07 fm. 

A main ingredient for the study of the equation of state is a good knowledge of the beta func- 
tions along the line of constant physics. In a direct evaluation one would measure the values of 
amp,mT,/mp,m^-Jm^ and fit those as functions of the three lattice couplings (3,Kud,i^s- However, 
the quark mass beta functions, adKi/da, are much smaller than the gauge part and the errors in the 
matrix inversion are often too large for an accurate determination. The authors of [113j therefore di- 
rectly fit the lattice couplings /3, Kud, i^s as a function of the meson masses using a third order polynomial 
ansatz in arrip, {niT^/mp), {niji-Jm^). As an example, the result for the lattice gauge coupling is shown in 
Fig. [25] (left), with similar plots for the other lattice couplings [113j . These fits can then be numerically 
converted to the required beta functions. 

The resulting trace anomaly is depicted in Fig. [25] (right). The thick error bars give an estimate 
of the systematic error due to variations of setting the scale by different quantities. Akima splines are 
used for interpolation and the numerical integration of the trace anomaly to yield the pressure. Starting 
point for the temperature integration is A,- = 16, for which T is low enough for the trace anomaly to 
be statistically compatible with zero. The authors also find that the variation associated with using 
different interpolation methods for the integration is smaller than the error bars. This is an interesting 
observation since, in the fixed scale method, simulation points for different temperatures can only be 
varied in discrete steps. The current spacing of temperature values could be halved by inserting odd 
values for A^,-. In addition, the temperature resolution can of course be improved by simulating on finer 
lattices in the future. 
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Figure 25: Left: Global fit for /3 as a function of m^a. Square symbols show coupling parameters in the 
CP-PACS/JLQCD study. Solid lines show global fit results at each simulation point with corresponding 
mp/rriT, and nirj^Jm^. Right: Trace anomaly, energy density and pressure for QCD with clover improved 
Wilson fermions, 771^^ ~ 620 MeV. Thin and thick vertical bars represent statistic and systematic errors, 
respectively. The curves are drawn by the Akima spline interpolation. From |113j . 



There are also preliminary results for the equation of state based on the twisted mass formulation 
of Nf = 2 Wilson fermions. Using also tree- level Symanzik improvement, the relevant formula for the 
trace anomaly is 

+^ mU]i-).,, - {2a,I'-ll + 2..^) {tn.rS>).^, ) . (169) 

Note that k, is tuned to Kc{P) to ensure maximal twist and the quark mass is tuned with the twist 
parameter fi. 

Preliminary numerical results are based on simulations reported in [112] as well as so far unpublished 
data, running on 32'^ x A^^ lattices with N^. = 6, 8, 10, 12, where the standard integral method is used. 
The pion mass is m^r ~ 400 MeV. The required beta function is obtained from interpolations of the 
chirally extrapolated vq parameter provided by the T = simulations of the ETMC collaboration |116] . 
This is shown in Fig. [26] (left). Similarly, the vacuum expectation values needed for renormalisation are 
obtained from interpolations of ETMC data. The resulting trace anomaly is shown in Fig. [26] (right). 
While it is too early to draw conclusions at this stage, scaling appears to set in at Nt->8 which would 
permit a continuum limit once enough statistics has been collected. 
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9 Comparing results from different discretisations 

Once a controlled continuum extrapolation has been performed, all lattice actions must, for fixed physi- 
cal parameters and observables, give the same answers, lest a particular lattice formulation is fundamen- 
tally flawed. Working with different discretisation schemes therefore is a necessary and powerful check 
of the validity of the results. It is however natural to also compare results even before the continuum 
limit has been taken. Different actions have different cut-off effects. Notable discrepancies are therefore 
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Figure 26: Left: Interpolation of the chirally extrapolated scale parameter tq from the ETMC collabo- 
ration |116] . Right: Trace anomaly for Nf = 2 twisted mass fermions with ~ 400 MeV. 

an indicator that the quantity under investigation has not yet reached its continuum value in at least 
one of the studied schemes. 

9.1 Staggered vs. staggered fermions 

Comparisons of the trace anomaly between different staggered fermion actions were made in |103l 1117] . 
Note that in pure gauge theory, A^,- = 6,8 with an unimproved action proved to be fine enough for 
a continuum extrapolation in the ~ scaling region, passing later tests by improved actions, as 

discussed in Sec. |6l This might have been the reason for an overly optimistic approach to dynamical 
simulations. A comparison between p4, asqtad, stout and HISQ/tree actions from [117j is shown for 
Nr = 8 in Fig. |27] (left). Striking differences are apparent. Note that the pions have physical masses 
for the stout data, but are slightly heavier for the other discretisations, m^r ~ 160 MeV. However, 
the quark mass dependence quickly diminishes for T > T^, since there are no more hadrons and the 
screening masses in the same quantum number channels are dominated by ~ 27rT, compared to which 
the vacuum masses of the pions are small. Hence, the quark masses should only affect the data below Tc. 
By contrast and consistent with the calculations of Tc alluded to earlier, features like the rising flank and 
the peak position of the hotQCD results appear shifted towards higher temperatures. Much larger than 
this shift, however, is the difference in the peak height, which at the maximum corresponds to almost 
50%. If it is not the difference in quark mass, the difference must be accounted for by discretisation 
effects. It is interesting to note the evolution of the different discretisations with N^- and compare with 
the corresponding evolution of the stout data shown in Fig. However, at least for the A^^ = 6, 8 data 
shown, neither the HISQ/tree nore the asqtad data decrease the peak in the trace anomaly noticably. 
It will be necessary and highly interesting to see A^^ = 12 data with the HISQ/tree action to clarify 
this discrepancy. 

Similarly, differences are observed when looking at the chiral condensate and the strange quark num- 
ber susceptibility. This is shown in Fig. |27] (right) and Fig. |28] (left), respectively. In these observables a 
very clear trend can be seen. The data of both, the asqtad and the HISQ/tree action, move significantly 
towards the stout results with increasing A^,-, i.e. on finer lattices. A smaller discrepancy is observed 
in the renormalised Polyakov loop. Fig. [28] (right). Note the gradual decrease in cut-off effects in these 
three observables. The light quark condensate A^^^ is particularly sensitive to the chiral properties of 
the action. These are dominated by the pion sector and the associated cut-off effects by taste splitting. 
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Figure 27: Left: Trace anomaly computed by different actions on A^^ = 8. Data for p4fat and asqtad 
are from [171 1101] . HISQ/tree from |117J . those for stout correspond to the continuum estimate based 
on A^^ = 6,8,10 in |103j . From |117j . Right: The subtracted chiral condensate for the asqtad and 
HISQ/tree actions [63] with mi = ms/20 compared with the continuum extrapolated stout action 
results [103] . The temperature T is converted into physical units using ri. From [63] . 

Accordingly, the actions with the smallest taste splittings, i.e. the stout and HISQ/tree actions are 
expected to do best for this observable. By contrast, the Polyakov loop is a pure gauge observable and 
couples only indirectly to the sea quark effects. It is therefore less sensitive to taste splitting. For all 
observables, the scaling seems to be improved when the scale is set with Jk rather than ri [661 163] . 

Finally, it is instructive to compare the numerical results below with the hadron resonance 
gas prediction. In Sec. 13.121 we have seen that the hadron resonance gas is more than just a model 
description, but becomes increasingly accurate the lower the temperature, and hence the lattice gauge 
coupling at fixed N^-. We have also seen how to adapt the hadron resonance gas to finite lattice spacings, 
by including the taste sphttings as well as the hadron spectrum shifted by finite cut-off effects. Sec. 14.91 
Following [To], the BMW collaboration performed a comparative analysis of hadron resonance gas 
descriptions with and without inclusion of cut-off effects [68]. This is shown in Fig. [291 Up to T ~ Tc, 
the continuum extrapolated stout data are in complete agreement with the hadron resonance gas using 
the physical continuum spectrum. On the other hand, similar agreement is achieved for the asqtad 
action when appropriate cut-off effects for the hadron masses are included in the hadron resonance gas 
description, corresponding to the "distorted" curves in Fig. [291 

Let us attempt to draw some conclusions. If we take agreement with the hadron resonance gas for 
a criterion, we obtain indeed one explanation for all discussed observables. The trend for the p4 and 
asqtad actions to push features of the curves to slightly larger temperatures is then due to cut-off effects 
manifested in the hadron spectrum and mainly in the pion splitting. On the other hand, the cut-off 
effects in the dispersion relation appear to matter less in the transition region. This is not surprising, 
their improvement effect is relevant for the ideal or weakly interacting gas and hence at very high 
temperatures. The fact that the equation of state up to almost Tc is described by the hadron resonance 
gas shows to the contrary that it is strong coupling physics that is at work here, and hence the cut-off 
effects on the hadron spectrum that matter most. For this reason the stout and HISQ/tree actions 
perform best in this range of temperatures. (Note that an improvement scheme that is complete to 
some order in the lattice spacing works of course on all scales). All of this is of course also consistent 
with the notion of a strongly coupled quark gluon plasma just above T^. Finally the observation that, 
for all observables considered, the stout data for > 8 appear to sit on top of each other within 
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Figure 28: Left: Strange quark number susceptibility for ^ 160 MeV with the asqtad and HISQ/tree 
actions compared to the continuum extrapolated stout results at physical masses |103] . Right: Renor- 
malized Polyakov loop for the asqtad and HISQ/tree actions compared with the continuum extrap- 
olated stout result |1U3] . From [53j . 



small errors indicates that, for the temperature range considered, they have reached the scaling regime 
and a reliable continuum extrapolation can be performed once the Nr = 12 data are complete. This 
constitutes the last step in the determination of the equation of state within that discretisation scheme. 
It is now most desirable to confirm this with an independent discretisation, and the HISQ/tree action 
is the most promising candidate. 



9.2 Staggered vs. Wilson fermions 

Wilson and staggered quarks have different systematic errors and hence make for a valuable one-to-one 
comparison of results in physical units. In the continuum limit, both must give the same results of 
course, if there are no fundamental problems with either discretisation. Conversely, deviations at finite 
lattice spacings give valuable insight into systematic errors and limitations, since in this case we are not 
only comparing cut-off effects of the same kind with different size, but really two conceptually different 
ways of dealing with the doubling problem. The difficulty is to find simulations performed at the same 
physical masses so as to be directly comparable. So far the equation of state itself is not yet available for 
this purpose, but the behaviour of the chiral condensate and the Polyakov loop through the transition 
have been studied by the BMW collaboration in a direct comparison between Wilson and staggered 
quarks |2SIEZ]. 

For the Wilson simulations they use the fixed scale approach, employing a tree-level Symanzik 
improved gauge action. Sec. 14.21 The Wilson fermion action has a tree-level clover improvement term, 
Csw = 1 Sec. I4.3[ and employs stout smeared links with six iterations and smearing parameter p = 0.11, 
Sec. 14.71 The scale was set by mn = 1672 MeV and the line of constant physics was determined by 
the mass ratios in Table |6l utilising data from previous spectrum calculations at zero temperature 
|118[I119[I120] . This translates to 171-,^ ~ 545 MeV and ttik ~ 614 MeV. On the other hand, the staggered 
results entering in the comparison are those obtained at fixed N^, following the same procedure to fix 
the lines of constant physics as in \V21\ EH EB] , cf . Sec. 17.41 Essentially the beta function for the 
strange quark was determined while keeping the bare quark mass ratio fixed at mi/mg = 2/3. The 
resulting ratios mT^/mQ,mK/'mQ were then found to be within 3% of the Wilson results. In order to 
renormalise the Polyakov loop, the authors of [26t [27] simply fix a value L* = 1.2 at a temperature 
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Figure 29: Left: Strange quark number susceptibility from continuum extrapolated stout data [68] at 
physical masses, compared to HISQ/tree and asqtad |17l 1117] . Also shown are hadron resonance gas 
predictions using physical and cut-off-distorted spectra. Right: Renormalized Polyakov loop for the 
asqtad and HISQ/tree actions [63] compared with the continuum extrapolated stout result |103] . From 



/3 






ms/rriud 


arriQ 


a [fm] 


3.30 


0.332(3) 


0.373(3) 


1.529(2) 


1.16(1) 


0.139(1) 


3.57 


0.319(6) 


0.359(4) 


1.531(2) 


0.777(9) 


0.093(1) 


3.70 


0.326(5) 


0.369(5) 


1.531(3) 


0.586(8) 


0.070(1) 


3.85 


0.314(7) 


0.358(6) 


1.528(4) 


0.480(8) 


0.057(1) 



Table 6: Spectroscopy from zero temperature Wilson simulations. The quark mass ratios are nis/mud = 
[Irn]^ — ml) /ml. The lattice spacings are set by = 1672 MeV. 

T* = 0.143mn. The renormalised Polyakov loop is then 

= ilk)) ' ■ 

A comparison between results from staggered and Wilson fermions is shown in Fig. |30] [26]. For the 
Polyakov loop good agreement is observed in the entire range. For the chiral condensate, the Wilson 
data corresponding to the two larger /3-values nicely fall on the curve of the staggered data, while the 
two coarser lattice spacings show large deviations that get stronger with temperature, as expected in the 
fixed scale approach. Note that the Polyakov loop appears to be less sensitive to discretisation effects 
than the chiral condensate, which is presumably because the latter is sensitive to the chiral symmetry 
breaking of the Wilson action at finite lattice spacing. 

Interestingly, this cut-off effect of the Wilson action can be largely reduced by employing the defini- 
tion Eq. ( 159|) for the renormalised chiral condensate. For an (9(a)-improved action, its corrections start 
only at O(a^), with leading order corrections cancelling between m/j and {ipip)^. Of course, this defini- 
tion requires in addition the definition and calculation of Za and mpcAc- With the combined tree-level 
clover and stout link improvement this appears to be almost achieved and a much improved scaling 
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Figure 30: Comparison of the renormalised Polyakov loop (left) and chiral condensate (right) between 
simulations using staggered and Wilson fermions. From [26]. 
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Figure 31: Continuum extrapolated results for the renormalised Polyakov loop (left) and the chiral 
condensate (right). Perfect agreement is found between Wilson and staggered fermion discretisations. 
From 1271. 



behaviour for the chiral condensate of the Wilson data is observed [27]. This allows for a continuum 
extrapolation of both types of discretisations over the whole temperature range, providing perfectly 
agreeing continuum limits, Fig. [3TJ 

We may thus conclude that, for pions with > 500 MeV, thermodynamics simulations using the 
staggered or Wilson discretisation schemes indeed arrive at the same continuum limit, i.e. the different 
lattice artifacts disappear from each formulation as desired. However, there is one caveat left prohibiting 
us from drawing the same conlusion for all quark masses. Concerns about a potentially wrong continuum 
limit of staggered fermions employing the rooting trick are to a large extent based on the U{1)a anomaly 
and its relation to the fermionic zero modes [21]. This potential problem is not an issue for the relatively 
heavy quarks considered here, but becomes relevant when quarks are light such as in the physical case. 
It is thus mandatory to continue such comparisons until Wilson simulations either reach or can be 
extrapolated to the physical mass values. 
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10 The equation of state at finite density 



Heavy ion collisions proceed at a given baryon number, depending on the choice of nuclei. Within the 
grand canonical potential, this corresponds to a non-vanishing chemical potential for baryon number, 
which can also be tuned by the collision energy of the ions. A long term goal for future facilities like 
FAIR at GSI (Darmstadt, Germany), is to also probe cold high density matter that is relevant for the 
interior of compact stars. Unfortunately, few first-principles predictions exist for these regimes, because 
of the so-called sign problem of lattice QCD. For a complex quark chemical potential, the determinant 
satisfies 



which is complex unless Re /i = 0. However, a complex fermion determinant cannot be interpreted as 
a probability measure for importance sampling and the evaluation of the path integral by Monte Carlo 
methods is spoiled. On the other hand, simulations at imaginary fi or finite isospin chemical potential, 
/^M = — /^d, feature a positive determinant and can be simulated in the same way as at /i = 0. 

It should be stressed that this is merely an algorithmic problem, the partition function and physics 
remain completely well defined. The same problem appears also in condensed matter physics, whenever 
fermionic systems with a net density of particles vs. anti-particles or particles vs. holes are considered. 
Existing solutions to sign problems with cluster algorithms |122] . world line representations |123] or 
complex Langevin simulations |124j could not yet be adapted to full dynamical QCD. Available methods 
for QCD are reweighting techniques [125^, Taylor expansions about fi/T = |126t 1127] . simulations at 
imaginary chemical potentials followed by analytic continuation [128^ 1129] and use of the canonical 
ensemble |130j . All of these introduce (partly related) additional systematic errors and are valid only 
for /i/T < 1. For a pedagogical introduction, see [3l]. Because of technical difficulties and computational 
effort associated with these calculations, most are restricted to comparatively coarse and small lattices. 
Recent reviews are given in [HI I13H I132[ 1133] . Here we merely collect some results for the effect of finite 
baryon density on the equation of state. 

First qualitative results for the dependence of the equation of state on quark chemical potential were 
produced using reweighting techniques in |134] . for coarse lattices with A^,- = 4 employing standard 
staggered fermions corresponding to m,r ~ 400 MeV. Similarly, the Taylor expansion approach with 
cumulant approximations to the fermion determinant was tested on A,- = 4 lattices with two fiavours of 
Wilson fermions with m^r ~ 500 MeV [135] . In the meantime, there are more realistic calculations, all 
of them based on the Taylor expansion approach which we thus review here, after some more general 
considerations. 

10.1 Qualitative results from Taylor expansion and imaginary fi 

The study of the finite density modifications of the equation of state based on a Taylor expansion is 
straightforward. Essentially, corrections to the zero density equation of state are calculated order by 
order in sufficiently small chemical potentials. Let us begin with some qualitative considerations for the 
theory with Nf = 2 [136]. The pressure is written as a Taylor series 



where fnm are the Taylor expansion coefficients, which vanish for odd (n + m) because of CP-symmetry, 
and the zero density contribution has been subtracted. When the quark masses are equal, one has 
fnm = fmn- The Standard procedure is to calculate the coefficients at /x = 0. However, in [136] it 
was remarked that the /nm's are also related to derivatives Xij of the pressure measured at non-zero 
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n,m=l 



60 




Figure 32: Left, Middle: Taylor coefficients fn and /2o- Right: x^/dof of various polynomials as a 
function of T /T^ with fitting range a/i/ is 0.0-0.24. From |136] . 



chemical potential by 

n=t,m=j 

At zero density Xij = fijT'^~^~\ while at non-zero densities Xij includes higher order fnm terms, and does 
not vanish for odd {i + j)- The proposal of [136j is to evaluate the Xij imaginary or isospin chemical 
potential, where simulations can be carried out, by fitting all of them simultaneously to the polynomial 
expansion Eq. (I173p . In practice, imaginary chemical potential /i = i^i was used and eight different 
Xij evaluated up to i -|- j = 4. The data were then ffited to the corresponding polynomial expansions 
Eq. f ll73p truncated to a given order {n + m), which gives the corresponding fnm- Knowledge of the 
coefficients fnm then allows to reconstruct the thermodynamic potentials. 

In the confined phase the equation of state is well described by the hadron resonance gas. In this 
case the expansion in chemical potentials reads [137] 

= G[cosh(^)-l]+i?[cosh(^)cosh(^)-l] (174) 



+ iy[cosh(^) (^cosh(^) + cosh(:^)J - 2], 

where G, R and W are constants related to the hadron spectrum which can be translated again to 
the fnm- Quark and isospin chemical potentials fig and fijs are defined as fig = {fiu + /^d)/2 and 
/i/s = - lJ'd)/'2. respectively. 

To test those propositions, simulations have been performed for Nf = 2 flavours of standard stag- 
gered quarks on exploratory 8^ x 4 lattices, with a bare quark mass am = 0.05, with a range of imaginary 
chemical potentials a/Zj = 0.0, . . . 0.24 and a temperature range T/T^ = 0.83, . . . 2.0 |136j . Fig. 13^ shows 
examples for two Taylor coefficients extracted by fitting a sixth order polynomial. They are found to 
agree well with coefficients computed directly as derivatives, but the results from ffited polynomials 
have smaller errors. The main finding of the investigation concerns the possibility to capture the physics 
accurately with truncated polynomials. In Fig. [32] (right) the x^/dof of the ffis is shown as a function 
of temperature. A fourth order polynomial does not allow for good ffis anywhere in the temperature 
range, and even a sixth order polynomial has difficulties around Tc. The conclusion is that for a reliable 
description one needs eighth order polynomials at least, which is the best currently available information 
on coarse lattices within the Taylor series approach [138] . 
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10.2 Equation of state to order fi with almost physical quarks 

The Taylor expansion approach is followed by the MILC and hotQCD collaborations for Nf = 2 + 1 
flavours, with chemical potentials for the light and strange quarks, fii^s/T, respectively |107] . The 
pressure and trace anomaly are now sums in these variables. 



n,m=0 



E ^^rn{T) ^ ^ , (175) 



rp4 

" n,m 



The bar over the chemical potential indicates that they are measured in units of MeV rather than in 
units of temperature. The partition function is an even function of chemical potential due to its CP 
symmetry, i.e. n + m is even. The coefficients Cnm, bnm are the appropriate derivatives evaluated at zero 
density. Since all /^-dependence resides in the determinants, what is required are expressions like 
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Tr M-^-—M-^ 



(177) 



These are traces of local operators which can be evaluated using random noise vectors. Explicit expres- 
sions for the coefficients can be found in |139] . 

The coefficients have been computed through /i^ using the configurations generated with the asqtad 
staggered fermion action for physical strange quark mass and m^/m^ = 0.1, ~ 220 MeV, i.e. those 
entering the zero density calculations covered in Sec. 17.31 Examples for some of the coefficients are 
shown in Fig. |33] (top) for lattices with A^^ = 4, 6 |107] . Given previous experience at zero density, it 
is not surprising that in some coefficients there are differences between the Nr = 4: to Nr = Q results, 
indicating the sizeable discretisation errors present at N^- = 4. 

Once the coefficients are at hand, the change AX = X{fii^s) — X(0) in the thermodynamic potentials 
can be computed. Note that since the c„i(T) terms are non-zero, a non- vanishing strange quark density 
Us is induced even for /i^ = 0. Conversely, in order to have zero strange quark density, one needs to tune 
the strange quark chemical potential /is(T, /i/). The corresponding results for the pressure are shown in 
Fig. [33] (bottom left). As in the coefficients, cut-off effects are clearly visible in the equation of state. 

In heavy ion collisions, once the plasma has thermalised, it expands to cool at constant entropy. To 
get closer to the experimental situation, it is then interesting to consider the isentropic equation of state, 
for which the ratio of entropy to baryon number is held fixed. Isentropes in heavy ion collisions are 
trajectories in (/i/, /i^, T)-space with Ug = 0, which are characterised by different values for s/hb = C. 
For AGS, SPS, and RHIC, s/ub ~ 30,45, and 300, respectively |140j . For a two flavour theory, they 
were flrst calculated on the lattice in |141] . The pressure evaluated along those isentropes is shown in 
Fig. |33] (bottom right). Interestingly, the discretisation effects appear smaller in these quantities. The 
authors explain that this is due to the fact that the zero density contribution, i.e. the flrst term in the 
Taylor series, is largest for these quantities and there the difference between Nr = 4,6 was observed to 
be small [98] . However, it should be kept in mind that these lattices are too coarse for scaling and this 
can not be taken as an indication of the proximity to the continuum, as we have seen in Sec. |9l 
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Figure 33: Top row: Some coefficients Cnm in the Taylor expansion of the pressure as a function of 
temperature. Arrows indicate the respective (massless) Stefan-Boltzmann hmits. Bottom row: Change 
in the pressure due to the non-zero chemical potentials (left) and pressure along different isentropes 
(right). From |107] . 

10.3 Equation of state to order fi^ for physical quark masses 

Recently, the finite density corrections to the equation of state have also been evaluated by the BMW 
collaboration |142j . making use of their /i = configurations generated with the stout link improved 
staggered action at physical quark masses for various lattice spacings. The results discussed here are 
the finite density extension of those covered in Sec. 17.41 Specifically, they look at the leading order 
correction, which is the second derivative with respect to chemical potentials. 



_ TJ_ dHogZ 
VT^ dfiidfij 

The trace anomaly to quadratic order is then 



(178) 



Chemical potentials were fixed such as to keep strange quark number density to zero, Ug = 0, 

= IJ'u = fJ-d, jJ's = -2^fiu ■ (180) 

X2 

Alternatively, one may assign the same chemical potential for baryon number to all quark flavours. 

The required derivatives were calculated in the temperature range 125 MeV < T < 400 MeV. 
We recall that configurations on A^^ = 6, 8, 10, 12 were available for this purpose, which have been 
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Figure 34: Left: Difference between tlie pressure at /i > and /i = 0. Right: Trace anomaly for 
non-zero /il- From |142] . 



supplemented by = 16 [143] . Hence all intermediate quantities can be extrapolated to the continuum. 
For lower temperatures the partition function for the hadron resonance gas was used in order to obtain 
the appropriate derivatives. The following figures are therefore to be interpreted as continuum results 
at the physical point. 

Fig. EH (left) shows the difference in the pressure due to finite chemical potential for the light 
quarks. Again perfect agreement with the continuum hadron resonance gas is observed roughly up to 
the transition temperature. The change in the trace anomaly due to chemical potential fii = 400 MeV 
is shown in Fig. |31] (right). In a similar fashion, results for the pressure, the energy and densities as 
well as the speed of sound are compared for /i = 0, 400 MeV in |142j . 

In a second step, the isentropic trajectories corresponding to constant entropy over quark number, 
S/N, are located. Fig. [35] (left). Note that there are significant deviations from leading order pertur- 
bative predictions. Furthermore, in the low temperature region T < 130 MeV, the isentropes require 
chemical potentials too large for the Taylor expansion to be valid, so these results are to be taken with 
care. The thermodynamic functions evaluated along the isentropes are shown in Figs. [351 EEl As in the 
case of zero density, the authors of [142J provide convenient parametrisations that reproduce the data 
over the whole range < T < 400 MeV, 



X2{T) 



/o-[tanh(/i-t + /2) + l] 



■/3-[tanh(/4-t + /5) + l], 



:i811 



with the appropriately fitted constants given in Table [3 The trace anomaly for small chemical potentials 
follows then from Eq. ( I179p . and the pressure by T- integration. 
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Table 7: Parameters for the functions fll81|) . 



Finally, a comment on the reliability is in order. The results discussed in this section are valid 
in the continuum and for physical quark masses, but only up to order /x^. In order to check the 
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Figure 35: Left: Trajectories of constant S/N^ on the phase diagram. The hnes in the upper part of the 
plot represent the leading order perturbative expressions. Right: Trace anomaly along the isentropes. 
From [Wj. 
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Figure 36: Pressure (left) and energy density (right) along the isentropes. From [142] . 



systematics, also the subleading coefficient X4 was checked in |142] . The statement there is that the 
{Hb/TY contribution to the pressure makes up < 10% of the leading order contribution for fis/T < 2 
and <35% up to Hb/T < 3. Hence, good accuracy is only achieved at high temperatures around and 
beyond the transition. We should also keep in mind the finding of Sec. IIU.H according to which a 
good description of the functional /z-behaviour takes several more but the leading coefficient. Future 
calculations of higher order contributions are thus necessary before the //-dependence can be regarded 
as controlled. 



11 Summary and Outlook 

After many years of intense studies and improvements on the theoretical, computational and algorithmic 
level, we are converging towards a fully non-perturbative and quantitative understanding of the QCD 
equation of state at finite temperatures and zero density. Lattice studies using the staggered fermions 
with stout link improvement give accurate predictions for physical masses and temperatures 125 < T < 
1000 MeV, a complete continuum extrapolation based on four lattice spacings in this temperature range 
is around the corner. We thus have a fully non-perturbative prediction for QCD in the continuum. An 
independent check of these results with a different staggered discretisation is also on the way. A last 
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1 1 
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Wilson 
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Wilson 
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Reference 


Section 
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220 MeV 
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156 MeV 
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trS 


stout staggered 


135 MeV 


0.083 fm 


T ^ 100 - 500 MeV 
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El 


RG 


Wilson clover 


620 MeV 


0.07 fm 




[113J 
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160 MeV 


0.125 fm 
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Gauge action 


Fermion action 




Smallest a 


charm 


Reference 


Section 


asqtad 


asqtad 


220 MeV 


0.17 fm 


quenched 


[107J 


[73] 


trS 


stout staggered 


135 MeV 


0.1 fm 


dynamical 


[105] 


1731 








Nf = 2 + 1,12^0 








Gauge action 


Fermion action 




Smallest a 


order in /x 


Reference 


Section 


asqtad 


asqtad 


220 MeV 


0.17 fm 




m 


110.21 


trS 


stout staggered 


135 MeV 


0.063 fm 


/i^, cont. limit 


m 


110.31 



Table 8: Some phenomenologically relevant simulations of the equation of state covered in this review. 
trS=tree-level Symanzik, tad=tadpole, RG=fixed point or Iwasaki, p4, asqtad all refer to improvement 
schemes explained in Sec. HI The smallest lattice spacing a is measured at T 200 MeV. All tabulated 
Nf = 2 + 1 simulations work with the physical strange quark mass. 



source of potential systematic error concerns the question whether rooted staggered fermions have a 
fundamental flaw in their continuum approach. This will be answered by the corresponding Wilson 
simulations as discussed in this review. Extending those comparisons down to the physical light quark 
masses will require a few more years. Furthermore, people have started thermodynamical simulations 
using dynamical overlap fermions which have the correct (lattice) chiral symmetry |1441 1145] . However, 
it will take much longer for those to reach the same level of maturity. 

At low temperatures the lattice results for the equation of state match smoothly to the hadron 
resonance gas predictions, which in itself is a non-trivial finding supported by first principles calculations 
based on strong coupling expansions. Together with the fact that for temperatures up to 1 GeV 
one observes significant deviations ~ 10% from the ideal gas predictions, this is testament to the 
strongly coupled nature of the plasma around Tc, with significant contributions from strongly interacting 
soft modes up to much higher temperatures. Within pure gauge theory contact to the perturbative 
regime could be established and there is no fundamental problem extending this to the dynamical case. 
Therefore, we may conclude that a complete description of the continuum equation of state at zero 
density and all temperatures is within reach. 

Finally, some interesting extensions have also been achieved, though not yet at the same level of 
refinement. We have learned that the inclusion of the charm quark has a significant effect at high 
temperatures relevant for the LHC. A continuum limit for these predictions is also in reach within a 
few years. Moreover, significant progress has been made in the treatment of small baryon densities. 
Leading order /x^ modifications on the equation of state have been computed, and the technologies to 
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extend these to higher orders using Taylor expansion and simulations at imaginary chemical potentials 
or employing reweighting are available. A quick reference to all these developments is provided in Table 
[HI We may thus conclude that the quark hadron transition at zero and low baryon density, together with 
the associated equation of state, constitute the first fully non-perturbative, first principles predictions 
from the lattice for physical QCD under extreme conditions. 
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